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1.0  INTRODUCTION 

The  scientific  community  has  exploited  the  unique  capabilities  of  semi¬ 
conductors  to  produce  various  kinds  of  electrical  and  optical  components. 
However,  the  increasing  energy  levels  these  components  are  required  to  handle 
are  quickly  surpassing  the  materials1  stable  energy  handling  capabilities. 
Research  continues  to  expand  the  known  stable  energy  handling  modes  of  semi¬ 
conductors.  Still,  the  question  arises:  At  what  energy  level  or  under  what 
conditions  do  unstable  energy  handling  modes  occur  for  semiconductors? 

Yet,  this  question  reveals  a  deeper  dilemma  upon  examination.  Numerous 
efforts  have  approached  this  problem  utilizing  very  simple  experiments  and 
a  priori  failure  concepts  to  obtain  an  understanding  of  failure  phenomena 
associated  with  unstable  energy  handling  modes  that  can  be  applied  to  more 
complex  situations.  In  general,  this  type  of  approach  has  been  thwarted  by 
the  appearance  of  failure  phenomena  (Ref.  1)  that  are  not  consistent  with 
simple  theories  (Ref.  1).  It  Is  apparent  that  such  results,  from  simple 
experiments,  point  to  a  more  complex  mechanism  that  may  be  dlscemable  under 
more  sophisticated  modeling  techniques. 

Ensemble  Monte  Carlo  simulation  of  semiconductor  carrier  transport  has 
consistently  Improved  over  the  last  decade  to  the  point  of  achieving  excellent 
agreement  with  experiments.  Thus,  this  state-of-the-art  tool  Is  the  vehicle 
upon  which  this  effort  will  work  to  achieve  Its  objective:  The  modeling  of 
carrier  transport  In  semiconductors  under  high  optical  or  electric  field 
conditions  to  determine  the  physical  mechanisms  Involved  in,  or  triggering, 
unstable  modes  of  carrier  transport. 

To  develop  this  capability  at  the  Air  Force  Weapons  Laboratory  (AFWL) 
without  outside  help  would  be  a  major  project  requiring  a  large  number  of 
man-hours  and  computer  time.  Therefore,  the  decision  was  made  to  acquire 
the  basic  Monte  Carlo  program  and  upgrade  the  code  from  there. 

The  primary  source  for  the  computer  code  was  from  North  Carolina  State 
University.  This  program,  authored  by  Dr  M.A.  Littlejohn  and  Dr  J.  Hauser, 

Is  the  basis  used  by  all  state-of-the-art  semiconductor  Monte  Carlo  simulations 
and  is  the  basis  upon  which  this  effort  will  build. 
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Simulation  results  were  obtained  via  the  ensemble  Monte  Carlo  program 
operated  by  Or  Karl  Hess  at  the  University  of  Illinois  at  Urbana,  Illinois. 
This  program  is  the  most  advanced  simulation  of  its  type  in  the  world  and  pro¬ 
vides  check  cases  for  this  studies  code  development  and  modeling  effort. 

Dr  Howard  Volkin  was  responsible  for  this  efforts  code  upgrading  and 
running. 

Interpretation  of  program  results  were  directly  aided  by  Lt  Mark  Snyder, 
AFVIL/NTAT,  Dr  Erick  Kunhardt  of  Brooklyn  Polytechnical  Institute  and  Dr 
William  Portnoy  of  Texas  Tech  University. 
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2.0  BACKGROUND  OF  THE  PROBLEM 


This  study  was  initiated  by  the  results  enumerated  in  References  1,  2 
and  3  as  applied  to  finding  the  physical  mechanisms  involved  in  semiconductor 
device  failure  due  to  an  electrical  overstress.  These  studies  found  the  col¬ 
lective  behavior,  known  as  second  breakdown,  that  occurs  in  semiconductor 
devices  exposed  to  an  electrical  overstress  was  the  result  of  the  formation 
of  an  unstable  solid  state  plasma.  Formation  of  the  plasma  occurred  due  to  the 
transition  from  single-carrier  to  double-carrier  current  injection  or  from 
the  excitation  of  carriers  into  the  upper  conduction  bands  forming  a  single¬ 
carrier  plasma  dominated  by  a  particular  conduction  minima.  In  principle, 
the  physical  mechanisms  associated  with  plasma  formation  were  found  to  have  a 
strong  dependence  on  the  specific  material  band  structure,  and  device  archi¬ 
tecture.  Modeling  a  semiconductor  material  (without  a  device  structure)  for 
single-  to  double-carrier  injection  could  be  performed  for  specific  cases 
analytically,  but  higher  order  band  structure  effects  can  only  be  dealt  with 
on  a  fully  quantum  mechanical  level.  Thus,  a  Monte  Carlo  simulation  utilizing 
a  quantum  mechanical  representation  of  the  crystalline  lattice  (called  the 
pseudopotential  method)  was  seen  as  the  only  alternative  to  pursue. 

Additional  experimental  evidence  for  plasma  formation  was  found  on  exami¬ 
nation  of  laser-induced  electric  breakdown  in  solids.  Bloemberger  (Ref.  4) 
states  that  laser-induced  electron  avalanche  is  usually  the  mechanism  that 
determines  the  failure  threshold  in  pure  solids.  Multi  photon  ionization  of 
electrons  leads  to  the  formation  of  a  cloud  of  electrons  that  become  opaque 
to  the  laser  and  thus  highly  absorptive,  leading  to  thermal  failure. 

This  physical  process  is  reported  by  Bloemberger  to  be  responsible  for 
failure  at  the  surface  as  well  as  in  the  bulk  of  a  material.  Yet,  this  may 
be  true  only  if  physical  and  chemical  defects  are  not  prevalent  in  the 
material . 

The  concept  of  failure  put  forth  for  laser-induced  electric  breakdown  is 
similar  to  the  concept  presently  used  to  explain  electrically  stressed  failure 
of  materials.  A  conmon  scenario  would  be  as  follows:  A  material  Is  stressed 
by  an  electric  field  (optically  or  electrically)  Initiating  carrier  production 
and  movement.  If  the  field  Is  Increased  in  magnitude,  Impact  ionization  will 
occur  forming  double-carrier  Injection  (also  known  as  avalanche  injection). 
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As  the  number  of  carriers  increase  with  further  impact  ionization,  a  point  is 
reached  where  the  kinetic  energy  of  the  carriers  is  efficiently  imparted  to 
the  material  lattice  resulting  in  joule  heating.  The  joule  heating  increases 
the  kinetic  energy  and  number  of  carriers  involved  in  impact  ionization  form¬ 
ing  a  positive  feedback  cycle  (or  runaway  effect).  If  not  checked,  this  feed¬ 
back  cycle  increases  with  time  until  thermal  shock  occurs. 

Closer  examination  of  experimental  work  in  device  failure  (Ref.  1)  shows 
direct  evidence  of  fine  structure  in  the  carrier  cloud  formed  from  intense 
avalanche  injection.  This  evidence  shows  isotropic  device  current  collapsing 
to  filamentary  structures  which  grow  and  destroy  a  device  by  thermal  shock. 
Current-voltage  characteristics  associated  with  failure  also  exhibit  a  fine 
structure  that  is  not  descriptive  of  currently  known  device  physics.  Thus, 
this  explanation  of  material  failure  (avalanche  leading  to  runaway)  is  not  an 
entirely  accurate  representation  of  the  failure  phenomena  observed. 

Various  theories  have  attempted  to  explain  the  fine  structure  observed 
but  most  have  failed  to  explain  the  general  aspects  under  which  failure  occurs; 
i.e. ,  formation  of  current  filaments. 

This  study  will  attempt  to  simulate  the  general  aspects  occurring  when  a 
semiconductor  material  is  exposed  to  high  amplitude  electromagnetic  fields 
from  an  electrical  or  optical  overstress.  Particular  attention  will  be  paid 
to  those  aspects  that  may  lead  to  the  formation  of  a  solid  state  plasma. 
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3.0  THE  SOLID  STATE  PLASMA 


3.1  DEFINITION  AND  CRITERIA 

A  plasma  is  generally  defined  as  an  ionized  gas  of  electronically  quasi¬ 
neutral  particles  that  exhibit  a  collective  behavior  in  reaction  to  the  sur¬ 
rounding  environment  (Ref.  5).  The  term  quasi-neutral  refers  to  a  neutrality 
that  occurs  due  to  the  kinetic  action  of  the  particles  forming  an  electro¬ 
static  shielding  effect  between  like  charges.  This  is  known  as  Debye  shield¬ 
ing  and  forms  a  critical  concept  necessary  to  the  definition  of  a  plasma. 

The  term  "collective  behavior"  describes  how  a  plasma  reacts  as  a  single 
entity  to  the  surrounding  environment.  Individual  carrier  movement  in  a  plasma 
is  not  readily  discernable  at  the  microscopic  level  due  to  Debye  shielding, 
and  plasma  reaction  to  an  outside  force  Is  dominated  by  the  need  to  sustain 
equilibrium  and  keep  the  plasma  free  of  electric  potentials  or  fields.  The 
only  way  to  maintain  these  requirements  Is  by  reacting  as  a  single  unit. 

Criteria  have  been  developed  to  provide  proof  of  a  plasma  state  according 
to  the  above  definitions.  The  first  criterion  deals  with  the  Debye  length  or 
thickness  of  the  Debye  shielding  around  a  charged  carrier.  The  Debye  length 
is  defined  as 

(kT't  \  1 

J  <  tyaUm't  dimensions.  ( 1 ) 

where 

k  »  Boltzmann's  constant 
T  *  carrier  teaperature 
n  *  carrier  density 
e  *  electronic  ebarge 
<*  dielectric  constant 

Under  this  criterion,  the  Debye  length  must  be  much  smaller  than  the  dimen¬ 
sions  of  the  system  containing  the  plasma. 

This  criterion  gives  rise  to  a  second  criterion  based  on  the  premise 
that  If  the  Debye  length  is  much  smaller  than  the  system's  dimensions,  a  suf¬ 
ficient  number  of  carriers  should  exist  within  a  Debye  sphere  to  provide 
shielding,  or 

No  *  (2) 

3 
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where 

So  =  number  of  particles  within  a  Debye  sphere 
n  =  carrier  density 
XD-  Debye  length 

A  final  criterion  deals  with  the  plasma's  ability  to  possess  collective 
behavior  that  is  not  swamped  out  by  random  collisions.  To  fulfill  this  cri¬ 
terion,  the  frequency  of  an  oscillation  (o>)  in  the  plasma  times  the  mean  time 
between  collisions  (t)  must  be  greater  than  one,  or 

wr  > 1  (3) 

These  three  criteria  must  be  met  in  order  to  define  an  ionized  gas  of 
carriers  as  a  plasma.  Does  the  carrier  transport  in  a  semiconductor  possess 
the  qualities  necessary  to  fulfill  the  above  criteria  to  be  called  a  plasma? 

The  answer  is  yes,  as  demonstrated  in  the  following. 

A  semiconductor's  material  is  capable  of  transporting  a  charge  in  accord¬ 
ance  with  the  basic  principles  of  quantum  mechanics.  The  crystalline  potentials 
that  make  up  the  background  for  charge  transport  influence  carrier  motion  via 
quantum  mechanics  in  such  a  way  as  to  Impart  an  effective  mass  (Ref.  6)  that 
is  greatly  reduced  from  the  carrier’s  rest  mass  (Ref.  6). 

The  crystal  potential  Is  not  a  constant  —  it  is  a  tensor  and  changes 
with  orientation  of  the  crystalline  lattice  to  the  applied  field.  Therefore, 
the  effective  mass  is  a  tensor.  The  crystal  potential  also  varies  in  terms 
of  the  energy.  Thus,  a  carrier  at  2  eV  will  view  the  lattice  differently  than 
a  carrier  at  1.5  eV.  The  aspects  of  a  plasma  composed  of  2  eV  carriers  are 
different  than  those  of  a  plasma  of  1.5  eV  carriers.  But  how  are  these  plasma 
aspects  determined?  And  how  is  it  known  that  a  solid  state  plasma  exists? 

The  carrier  population  in  a  semiconductor  can  take  on  the  aspects  of  a 
plasma  due  to  the  Heisenberg  uncertainty  principle 

ApA (4) 


where 

A p  «  uncertainty  In  linear  nonentua. 

Ax  s  uncertainty  in  position, 
b  *  Planck's  constant  a  6.62  x  10~27  erg-sec. 


6 


y*  y  V •  V!  Vy* 
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The  temperature  of  carriers  in  a  solid  is  related  to  the  carriers' 
kinetic  energy  by 


or 


3  1  2 

-kT  —  -m  <  tr  > 

2  2 

_  .  m<  v*  ><  P2  > 

T‘— — iimr 


(5) 


where 

k  =  Boltzmann's  constant. 

T  s  carrier  temperature, 
m  =  carrier  mass. 

average  of  the  carrier  velocity  squared. 

As  the  number  of  carriers  increase,  Ax  decreases,  a p  increases,  and 
the  carriers'  temperatures  increase. 

The  intrinsic  carrier  concentration  in  a  semiconductor  will  give  a  rough 
estimate  of  the  average  distance  between  carriers  by  AiNn^.  Equation  4 
now  becomes 

(6) 

and  Equation  5  becomes 


* 

The  carrier  mass  (m  )  in  Equation  7  is  the  effective  mass  for  a  particular 
material.  Here,  this  quantity  is  treated  as  a  magnitude  as  long  as  the  major- 
ty  of  carrier  movement  is  along  one  crystal  orientation  only.  The  effective 
mass  for  most  semiconductors  Is  only  a  fraction  of  the  carrier's  rest  mass 
(Ref.  7).  Thus,  a  small  effective  mass’and  large  carrier  concentration  can 
create  a  large  carrier  temperature,  as  shown  in  Table  1. 


AFVIL-TR-85-116 


TABLE  1.  THE  UNITS  FOR  THE  EFFECTIVE  MASSES  (i.e.  m* 


mt  UNITS  ruR  THE  EFFECTIVE  MASSES  (i.e.  m*  ,  m  hh  and 

m*lh)  are  in  grams,  the  subscripts  stand  for  electron, 

HEAVY  HOLE,  AND  LIGHT  HOLE  RESPECTIVELY.  THE  TEMPERA¬ 
TURES  ARE  IN  KELVIN.  NO  APPLIED  FIELD. 


InSb 

1  X 

1017 

0.015 

0.39 

0.021 

10.6 

0.40 

7.5 

GaA a 

1  X 

10 

0.07 

0.68 

0.12 

2.3  x 

10-2 

2.3  x  10-3 

1.3  x  10 

1  X 

1018 

0.07 

0.68 

0.12 

10.5 

1.08 

6.14 

N-Ge 

1  X 

1015 

0.12 

0.34 

0.043 

6.1  x 

lO-2 

2.1  x  10-2 

0.17 

1  X 

10 18 

0.12 

0.34 

0.043 

6.14 

2.16 

17.1 

N-Si 

1  X 

10 1  ^ 

0.26 

0.52 

0.16 

6.1  x 

10*3 

3  x  lO-3 

1  x  10*2 

1  X 

10l8 

0.26 

0.52 

0.16 

2.83 

1.41 

4.6 

The  values  obtained  through  Equation  1  for  Debye  lengths  will  increase 
with  increased  carrier  concentration  and  carrier  temperature  and  possess 
values  as  follows  in  Table  2. 


TABLE  2.  CALCULATED  DEBYE  LENGTHS  FOR  THE  TEMPERATURES  DETAILED  IN 
TABLE  1. 


n(cnr3) 


Debye  Lengths  (in  ca) 
e  hh 


1 

X 

1017 

2.92 

X 

10-7 

5.7 

X 

10-8 

2.4 

X 

10-7 

1 

X 

10I4 

5.5 

X 

10-7 

1.2 

X 

10-7 

2.75 

X 

10-7 

1 

X 

1018 

7.7 

X 

10-8 

2.5 

X 

10-8 

2.7 

X 

10-7 

1 

X 

1015 

2.16 

X 

10-7 

1.2 

X 

10-7 

3.8 

X 

10-7 

1 

X 

1018 

6.8 

X 

10-8 

4 

X 

10-8 

1.2 

X 

10-7 

1 

X 

1014 

1.9 

X 

10-7 

1.3 

X 

10-7 

2.4 

X 

10-7 

1 

X 

10l8 

4.1 

X 

10-8 

2.8 

X 

10-8 

5.12 

X 

10-8 

These  values  can  be  inserted  in  Equation  2  to  determine  if  the  Debye 
shielding  is  effective  (Table  3). 

TABLE  3.  NUMBER  OF  CARRIERS  CONTAINED  IN  DEBYE  SPHERE  FOR  THE 
TEMPERATURES  DETAILED  IN  TABLE  1.  NO  APPLIED  FIELD. 

Number  of  carriers  within  a  Debye  Sphere  (cm-3) 

_ n(cm~3) _ |*de _ NPhh _ NDLh 

InSb  1  x  1017  1.0  x  10-2  7.8  x  10-5  5.7  x  10-3 

GaAs  1  x  1014  1.8  x  10"5  7.2  x  10"7  8.7  x  10"6 

1  x  1018  1.9  x  10-3  6.5  x  10-5  8.2  x  10-2 

N-Ge  1  x  1015  4.2  x  10-5  7.2  x  10“6  1.9  x  10"1* 

1  x  lO1®  1.3  x  10-3  2.6  x  10-4  7.2  x  10-3 

N-Si  1  x  10™  2.9  x  10"6  9.2  x  10“7  5.8  x  10*6 

1  x  10lB  2.8  x  10“^  9.2  x  10“5  5.6  x  lO™ 

These  tables  show  that  many  semiconductors  do  not  fulfill  all  of  the 
requirements  for  a  plasma  without  an  external  field  —  the  number  of  carriers 
within  a  Debye  radius  are  Insufficient  to  shield  a  carrier  from  external 
electromagnetic  forces. 

Equation  3  Is  the  final  criterion  to  determine  plasma  formation.  This 
equation  represents  the  quality  of  the  plasma,  or  sharpness  of  the  plasma 
oscillations.  If  the  plasma  oscillations  are  heavily  damped  by  collisions, 
the  plasma  Is  not  well  formed  and  reacts  sluggishly  to  external  forces.  The 
types  of  collisions  that  occur  In  a  semiconductor  must  be  examined  to  deter¬ 
mine  the  sharpness  of  a  solid  state  plasma. 

Collisions  In  a  semiconductor  fall  under  three  different  categories: 

(1)  collisions  with  the  lattice,  (2)  collisions  with  other  carriers,  and 
(3)  carrier  scattering  by  phonons  (lattice  vibrations).  Collisions  with  the 
lattice  are  of  prime  concern  when  the  lattice  captures  carriers  through  sta¬ 
tionary  traps  or  defects.  This  action  limits  the  carrier's  lifetime  or  ability 
to  contribute  to  carrier  flow  In  the  material.  Lifetime  is  dependent  on  trap 
density,  free-carrler  density  and  carrier  energy.  Collisions  with  other  car¬ 
riers  Is  not  Important  until  the  carrier  density  and  energy  are  extremely  large. 
Phonon  Interaction  is  only  Important  when  carrier  energy  Is  relatively  high 
compared  to  room  temperature.  Thus,  under  the  conditions  of  small  external 
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fields  or  temperature  gradients,  the  lifetime  gives  the  best  approximation  of 
he  time  between  collisions. 

The  ability  of  a  material  to  react  to  an  electromagnetic  wave  is  measured 
as  a  function  of  the  materials  dielectric  relaxation  time.  The  relaxation 
time  Is  the  time  needed  for  a  material  to  redistribute  carriers  to  maintain 
electrical  neutrality  (Ref.  8). 

The  dielectric  relaxation  time  also  determines  the  plasma  frequency  of 
the  material;  i.e.,  the  frequency  above  which  electromagnetic  waves  can  be 
transmitted  in  a  free-carrier  medium  according  to  dispersion  relations.  The 
relaxation  time^  is  found  from  the  following  equation 


'•  *  (sb) 


where 


£  a  relative  permltivity  times  the  free  space  permittivity 
a*  a  effective  carrier  mass 
n  a  free  carrier  density 

This  equation  yields  values  for  the  materials  listed  In  Table  4. 

TABLE  4.  PLASMA  FREQUENCY  CALCULATED  FROM  THE  DIELECTRIC  RELAXATION 
TIME  FOR  EACH  MATERIAL. 


Plasma  Frequency 
n  (cn*3)  (herts) 


1  x  101? 
i  x  101* 
1  x  10l8 
1  x  10*5 
1  x  1018 
1  x  10 141 
1  x  1018 


3.53  *  10 1 3 
6.03  x  10” 
6.03  x  10*3 
1.28  x  10^2 
4.07  x  101 3 
3.19  x  10” 
3.19  x  1013 


The  calculated  plasma  frequencies  show  that  plasma  reaction  (thus  sharp, 
ness  of  oscillation)  Is  much  faster  than  the  time  needed  ( *  1  ns)  to  capture 
carriers  at  traps  In  the  lattice.  Thus,  situations  can  arise  in  semiconduc¬ 
tors  where  most  of  the  criteria  for  definition  of  plasma  can  be  met  with  no 
injection  effects.  This  discussion  does  not  imply  that  solid  state  plasma 
exist  under  all  conditions.  However,  the  basic  criteria  should  be  rechecked 
under  different  physical  conditions  to  insure  compliance. 
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Application  of  an  electric  field,  optical  field,  or  external  heating  can 
raise  the  velocity  of  carriers  to  extremely  high  levels  (mix  10^  cm/s)  with 
a  correspondingly  high  temperature  («200  K) ,  a  larger  Debye  length  x  10‘5cm) , 
and  a  larger  N  (*5).  The  preceding  tables  and  the  figures  in  the  Appendix 
show  that  a  solid  state  plasma  can  satisfy  the  same  criteria  of  other  well  known 
plasma  forms. 

Applying  the  principles  put  forth  on  the  preceding  pages  to  more  specific 
cases  results  In  the  description  of  three  different  kinds  of  solid  state  plas¬ 
mas  depending  on  the  quantity  of  carriers  and  type  of  doping. 

3.2  COMPENSATED  PLASMAS 

A  compensated  solid  state  plasma  contains  equal  numbers  of  mobile  electrons 
and  holes.  This  type  of  plasma  is  likely  to  occur  under  Intrinsic  material 
conditions  at  low  level  Injection  regimes.  The  occurrence  of  avalanche  break¬ 
down  In  a  semiconductor  may  form  the  conditions  for  a  compensated  plasma,  if 
the  material's  band  structure  does  not  play  an  exceedingly  Important  role. 

This  type  of  plasma  can  be  weakly  or  strongly  Ionized  (a  large  or  small  carrier 
temperature) . 

3.3  UNCOMPENSATED  PLASMS 

This  plasma  forms  In  doped  semiconductors  where  the  concentration  of 
mobile  electrons  and  holes  are  unequal.  The  majority  of  mobile  charges  are  com¬ 
pensated  for  by  stationary  donor  or  acceptor  sites  In  the  crystalline  lattice. 

The  majority  carrier  density  under  these  circumstances  Is  relatively  high  and 
makes  formation  of  a  plasma  easy.  This  type  of  plasma  will  be  predominant  In 
semlconductdr  devices,  especially  under  high-level  Injection  regimes. 

3.4  UNCOMPENSATED  MULTICOMPONENT  PLASMS 

Once  free-carrier  energies  are  high  enough  to  reach  the  upper  conduction 
bands  of  semiconductors  In  large  nvaabers,  the  quantum  mechanical  aspects  of  the 
material  Influence  the  formation  of  a  plasma. 


Carriers  deposited  In  various  conduction  band  minima  will  possess  different 
dynamic  characteristics  than  other  carriers  in  other  minima.  Characteristics 
could  vary  by  substantial  amounts  in  terms  of  carrier  energy  and  effective 
mass.  These  differences  could  allow  plasmas  to  form  against  a  background  of 
other  mobile  carriers  that  do  not  share  the  same  energy  and  effective  mass. 
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This  type  of  plasma  would  be  highly  dependent  on  the  form  and  level  of  injection 
and  the  material's  specific  band  structure. 

3.5  NONEQUILIBRIUM  PLASMAS 

Situations  occur  In  semiconductor  materials  where  large  quantities  of 
charge  can  be  Injected  or  created  via  avalanche  breakdown.  This  charge  Is 
self  compensating  (equal  numbers  of  electrons  and  holes  are  created)  and  do 
not  require  compensation  by  Immobile  charge  centers.  This  type  of  plasma  Is 
often  strongly  Ionized  (due  to  high  electric  fields)  and  easily  populates  the 
upper  conduction  band  minima. 

The  four  types  of  solid  state  plasmas  outlined  In  this  section  would  also 
possess  different  forms  of  plasma  Instabilities.  The  Initially  formed  plasmas 
would  be  stable,  but  heating,  geometry  and  changes  in  external  electric  fields 
and  currents  would  force  the  plasma  Into  a  nonequilibrium  state. 

The  simplest  form  of  Instability  that  would  occur  Is  similar  to  the  plasma 
pinch  effect  Initiated  by  the  self-induced  magnetic  field  of  the  carriers. 
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4.0  EQUILIBRIUM  AND  INSTABILITIES  IN  PLASMAS 


Solid  state  plasmas  differ  from  gaseous  discharge  plasmas  in  many  aspects. 
In  a  simplistic  sense,  they  represent  extremes  of  the  definition  of  plasma. 
Differences  in  carrier  temperature  and  carrier  density  can  be  very  large  when 
comparing  thermonuclear  or  stellar  plasmas  to  solid  state  plasmas.  Yet,  all 
plasmas  share  the  concepts  of  equilibrium  and  stability  due  to  the  processes 
of  external  interaction  and  confinement. 


Since  the  semiconductor  material  is  the  confining  bottle  for  a  solid  state 
plasma,  problems  of  space-charge  buildup,  varying  current  flow,  high  electric 
fields,  and  nonlinear  carrier  dynamics  must  be  faced.  These  aspects  directly 
affect  the  plasma's  thermodynamic  equilibrium. 

If  the  plasma  Is  In  the  thermodynamic  equilibrium  with  Its  surroundings. 

It  Is  In  a  state  of  highest  entropy,  by  definition.  With  the  free  energy 
available  at  a  minimum,  the  plasma  remains  stable  unless  energy  Is  provided 
externally  to  excite  plasma  waves  (Increase  the  available  free  energy).  In  a 
semiconductor,  this  additional  free  energy  can  come  from  a  number  of  sources, 
but  Is  limited  In  general  to  heating  and  external  current  sources.  Plasma 
waves  use  extra  free  energy,  at  the  expense  of  the  plasma  carriers,  to  rapidly 
grow  and  dominate  plasma  behavior.  In  reaction,  the  plasma  assumes  motion 
to  decrease  the  available  free  energy  and  drive  the  plasma  back  toward  thermo¬ 
dynamic  equilibrium.  This  action  Is  known  as  an  Instability. 

Instabilities  are  classified  under  four  categories  depending  on  the  driv¬ 
ing  force  (or  free  energy)  Initiating  the  instability.  These  categories  follow 
from  basic  Instability  phenomena  found  In  gaseous  discharge  plasmas. 

4.1  STREAMING  INSTABILITIES 


When  a  current  is  driven  through  a  plasma,  various  species  of  carriers 
can  take  on  different  drift  velocities  with  respect  to  each  other.  Plasma 
waves  can  grow  by  Interacting  with  and  stealing  energy  at  the  expense  of  the 
carriers'  drift  velocity  resulting  In  Landau  damping. 

4.2  RAYLEIGH- TAYLOR  INSTABILITIES 


This  type  of  Instability  comes  from  forming  a  nonuniform  plasma  that  is 
driven  from  equilibrium  by  thermal  or  diffusive  forces.  A  comnon  example  In 
solid  state  plasmas  is  the  Gunn  effect. 
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4.3  KINEMATIC  INSTABILITIES 


When  the  velocity  distribution  of  carriers  in  a  plasma  are  not  Maxwellian, 
anistropies  can  develop  in  plasma  growth  oscillation  rates.  For  example, 
the  perturbation  of  a  plasma  changes  some  of  the  plasma's  potential  energy 
into  kinetic  energy  of  motion. 


4.4  UNIVERSAL  INSTABILITIES 

These  Instabilities  occur  due  to  the  way  a  plasma  is  confined.  If  a 
solid  state  plasma  expanded  very  rapidly  and  was  physically  confined  by  the 
material,  instabilities  would  occur  due  to  the  rapid  increase  in  free  energy 
available.  This  type  of  Instability  may  be  difficult  to  observe  because  of 
the  predominance  of  the  other  three  types  of  Instability. 

These  instabilities  normally  serve  to  dampen  oscillations  In  a  plasma, 
but  the  normal  driving  forces  (free  energy  sources)  in  semiconductors  act  in 
a  nonlinear  fashion  forcing  oscillation  growth.  To  control  these  instabilities, 
extraordinary  precautions  must  be  taken  to  control  a  semiconductor's  non¬ 
linear  dynamics.  This  often  Includes  a  liquid  nitrogen  bath,  quasi-steady 
state  electric  field  application,  and  application  of  an  external  confining 
magnetic  field. 

Without  precautions,  a  solid  state  plasma  can  rapidly  move  away  from 
thermal  equilibrium  with  the  lattice  with  the  Insulng  Instabilities  forcing 
material  damage  from  thermal  shock. 
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5.0  DEVICE  AND  MATERIAL  FAILURE 

The  postulation  of  the  formation  of  solid  state  plasma  in  semiconductors 
caused  by  an  overstress  provides  an  additional  clue  to  determining  semi¬ 
conductor  failure.  Since  ionized  solid  state  plasmas  move  away  from  thermo¬ 
dynamic  equilibrium  under  normal  device  operation,  the  determination  of  that 
departure  point  is  very  important. 

Previous  work  in  the  area  of  device  failure  attempted  to  model  failure 
in  a  very  simplistic  fashion.  Most  models  chose  to  relate  thermal  heating 
of  the  device  to  the  input  power  via  a  simple  empirical  relation.  These  models 
could  identify  gross  trends  in  device  failure  versus  input  pulse  length,  but 
sacrificed  the  accuracy  of  individual  device  failure  levels  in  favor  of  ensem¬ 
ble  trends.  Attempts  to  relate  this  Information  to  more  detailed  device 
physics  has  been,  to  date,  unsuccessful. 

Logically,  this  train  of  events  leads  researchers  in  device  failure  to 
examine  the  only  area  of  research  left:  The  detailed  physics  of  device 
failure. 

Second  breakdown  has  been  the  traditional  Indicator  of  ensuing  device 
failure  for  many  years.  The  phenomena  associated  with  second  breakdown  have 
also  been  well  documented  over  this  same  time  period.  The  evidence,  detailed 
in  many  reports,  points  to  the  formation  of  a  fine  structure  in  the  transport 
dynamics  of  the  materials,  that  Is  unprecedented  In  ordinary  knowledge  of 
carrier  transport.  Simply  stated,  the  Isotropic  current  before  second  break¬ 
down  suddenly  collapses  Into  a  high  current  density  channel  which  grows  and 
destroys  the  material  by  thermal  shock. 

The  current-voltage  characteristics  displayed  as  second  breakdown  occurs 
are  very  similar  to  curves  studied  for  the  transition  of  single  carrier  to 
double  carrier  Injection.  Thus,  It  becomes  apparent  that  a  strong  inference 
can  be  made  connecting  second  breakdown  to  solid  state  plasma  formation.  In 
particular.  It  appears  the  common  characteristic  of  current  filament  formation 
Is  similar  to  the  strongly  Ionized  plasma  instability  known  as  the  Bennett 
pinch.  This  Instability  Is  quickly  followed  by  rapid  heating  (also  a  conse¬ 
quence  of  the  Bennett  pinch)  that  could  easily  melt  a  semiconductor. 
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Other  forms  of  instability  can  develop  (as  noted  in  the  section  of  insta¬ 
bilities),  but  these  are  limited  primarily  to  specific  cases  of  device  archi¬ 
tecture,  overstress  characteristics  and  ambient  conditions.  The  simplest 
connection,  and  possibly  most  insightful  example  of  the  relationship  between 
second  breakdown  and  device  failure,  is  the  examination  of  the  solid  state 
plasma  under  a  Bennett  pinch.  This  example  will  be  used  to  show  that  the  key 
to  understanding  device  failure  and  second  breakdown  is  intimately  connected 
to  the  stability  and  characteristics  of  the  solid  state  plasma. 
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6.0  FI LAMENTATION  DUE  TO  PINCH 


6.1  THE  BENNETT  PINCH 


Since  the  primary  phenomenon  associated  with  second  breakdown  is  the 
formation  of  a  current  filament,  a  theory  of  second  breakdown  incorporating 
the  solid  state  plasma  must  explain  this  phenomenon. 

In  plasma  physics,  the  formation  of  a  filament  from  an  otherwise  isotropic 
plasma  is  called  the  Bennett  Pinch  (Ref.  9).  Basically,  the  flow  of  current 
carried  in  the  plasma  driven  by  an  external  electric  field  results  in  the 
formation  of  a  self-magnetic  field.  The  resulting  interaction  of  E  and  B 
through  the  Lorentz  equation  produces  a  radial  force  on  adjoining  plasma 
tubes  causing  contraction  (Fig.  1). 

Under  equilibrium  conditions,  the  plasma  carriers  attain  an  isotropic 
temperature  in  a  time  long  compared  to  changes  in  the  pinch  radius.  When  this 
is  not  true,  a  dynamic  pinch  must  be  used  to  model  the  effects. 

6.2  THE  EQUILIBRIUM  PINCH 

The  equilibrium  pinch  starts  with  the  general  fluid  equation  of  motion 
from  kinetic  theory. 


mn  (|j  ♦  (U.  v]u) 


»  qn  (T  ♦  Ux5)  -  v  •  7  ♦  7 


■Ul 


(9) 


where 


m  «  mass  of  carrier 
n  *  number  of  carriers 

U  ■  velocity 
q  »  charge 
E  *  electric  field 
B  *  magnetic  field 

P"  »  pressure 

P.  *  carrier  collision  factor 
l  r 


Under  equilibrium 


♦  (U  •  V)U  •  0 


so  that 


0  «  gnCF  *  17xB)  -  v  •  V  +  V. 
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Figure  1.  Once  a  sizable  current  Is  established 
the  J  x  B  force  will  cause  an  attrac¬ 
tion  of  flux  tubes  -  and  a  pinching  of 
current  lines. 
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Linearizing  7^  in  one  direction  so  that 
7 


results  in 


O  =  g«(  7  ■+•  V xBh~  VP  +  7. 


(12) 


Writing  Equations  11  and  12  out  for  electrons  and  holes 


0  =  -VP,  -  N,e(E  +•  o.xB)  +  P*h  (13) 

0  =  -VP*  -  udxfl)  f  Phr  (14) 

J  ~  S,evr  ■+  JVfcewh  ( 15 ) 

and 

(16) 

and  adding  Equations  13  and  14  gives 

Wp  =* IxB  .  (17) 


For  a  cylindrical  plasma  column,  the  change  of  pressure  with  radius  is 

- -J,(r)S,(r)  (18) 

as  shown  in  Fig.  2. 

Since  the  situation  Is  symmetric,  concern  is  primarily  with  one  direction  for 
current  flow,  so  that 

/.(«•)=  I'  J z[r')2**' dr  (19) 


Placing  Equation  20 


y 


Figure  2 


Cylindrical  coordinate  system  applied 
to  a  plasma  column  under  consideration. 
The  plasma  is  assumed  to  be  infinitely 
long  so  end  effects  can  be  ignored. 
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into  Equation  18  results  in 


4  TT'r 


dP(r) 

dr 


M../,(r)-/;(r)  =  - 


dr 


(21) 


Integrating  both  sides  from  r  =  0  to  R,  and  integrating  by  parts  on  the 
left-hand  side  (LHS)  of  Equation  21  results  in 


M../?(r) 


_  mo /;(/?) 
2 


dr  =  — 


(22) 


using  /,(o)  =  o.  next  gives 

4*VP(r)|*  -4*  f  2 xrP[r)dr  = 

/«  2 

P(r)  is  nonzero  and  finite  for  o <r<R  and  zero  for  R<r<oo. 
the  LHS  of  Equation  23  is  zero  and  P(r)  »  N{r)k(Tt  +TK)  gives 

-  4*  I  2 xr  MT,  +  T*)Af(r)  dr  = 

Jo  2 


(23) 


(24) 


let  N,  ■  N[r)2*rdr  be  the  number  of  carriers  per  unit  length  of  the  column, 
then 


/*(*)= 


(25) 


is  the  critical  current  (in  MKS)  required  to  initiate  pinch  for  a  specific 
temperature  and  carrier  number.  Equation  25  is  known  as  the  Bennett  relation. 
For  a  unit  bar  of  silicon  with  a  doping  level  of  1  X  1014  carriers/cn»3  /,  la 
about  0.87  mA.  Of  further  interest  is  the  distribution  of  PCr),  7,(r), 
and  B,(r )  from  the  Bennett  relation. 


To  Investigate  how  Equation  18  changes,  the  plasma  are  assumed  to  be 
ideal;  that  is,  they  possess  an  infinite  conductivity.  Thus,  the  current 
is  restricted  to  the  surface  of  the  plasma  cylinder.  From  the  law  of  magne¬ 
tostatics  Vxfl(r)  a  poJ{v).  SO 


J.{r)  =  —  j{rBv{r))=l-i-B'(r)+^- 

Mo'  dr'  no  dr  nof 


(2G) 
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which  placed  in  Equation  18  yields 


ipW  =  -*wiB,w-aw 


dr 


Mo 

1  d 
2 Mo^2 


Mor 


(27) 


(”2^(r)} 


Integration  .from  zero  to  r  arrives  at 

Since  P(r)  =  0  we  get 

which  when  replaced  in  Equation  28  yields 


(28) 


(29) 


with 


*.(»■)  *  fo  Ja^2wrdr  “d  fl*(r)  * 


(«•) 


(30) 


(31) 


P(r)  can  be  examined  if  J,(r)  is  known.  By  assuming  that  J,(r)  is 
concentrated  at  the  surface  (r  ■  R),  then  £*(r)  exists  outside  the  column  and 
P(r)  *  P(0),  so  pressure  Is  constant  (see  Fig.  3). 

Also,  since 


<P{,)  >m  7rpj„ 


(32) 

(33) 


combined  with  Equation  17,  the  first  term  on  the  right-hand  side  (RHS)  is 
zero.  Thus, 

<  PW  >.  (34) 
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ARBITRARY  UNITS 


p(r>7  f 

Jz(r) 

\ 

J 

(1  N^-B0(r) 

2  r/1 


Figure  3.  Distributions  for  P(r),  J  z(r),  and  B(r) 
for  a  current  sheath  at  t6e  edge  of 
the  plasma  column. 
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The  average  pressure  inside  the  column  (which  is  constant  for  this  case) 
is  balanced  by  the  magnetic  pressure  at  the  column  surface. 

If  a  solid  state  plasma  exhibited  an  infinite  conductivity,  the  column 
would  possess  ballistic  transport  qualities.  Although  it  is  possible  to  attain 
ballistic  trajectories  in  a  semiconductor,  the  necessary  conditions  are  quite 
severe.  Thus,  ballistic  transport  is  unlikely  to  dominate  current  transport 
processes  under  normal  conditions.  It  is  possible  that  the  carrier  concen¬ 
tration,  carrier  mobility  or  transit  times  could  reach  levels  where  the  con¬ 
ductivity  was  quite  high  compared  to  other  portions  of  the  device  (resembling 
a  strongly  ionized  plasma).  Current  could  then  be  monotonically  increasing 
with  distance  from  the  center  of  the  cylinder  (this  has  been  postulated  in 
Ref.  10  and  the  above  derivation  seems  to  support  the  idea). 

Another  possible  distribution  is  Jt(r)  *  constant  across  the  entire 
cylinder  and  zero  at  r  *  R.  In  this  case, 

»•<'>•!&  <«> 

Substitution  of  this  value  of  B#(r)  into  Equation  30  generates 


a  parabolic  dependence  of  pressure  on  radius  (Fig.  4).  This  example  would 
be  found  for  a  weakly  ionized  plasma. 
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To  get  a  better  understanding  of  J(r)s  possible  distribution,  refer  to 
Bennett's  original  study  (Ref.  9).  In  that  study,  Bennett  found  the  distribu¬ 
tion  or  N(r)  which  relates  to  J(r)  via 


J,(r)  ■  JV(r)«tv,  +  Vfcj 


(37) 


which  relates  to  P(r)  and  B(r)  by 

it 

P[r)  «  N(r)k{Tt+  Tk) 


Tk)  -  (0«(v.  *  •*)*•(”) 


Multiplying  by  r/N(r)  and  taking  the  derivative,  replacing  with 

£<•*.(')  -  *(r)«(v.  ♦  *k)  - 
results  in 


s{(sr)5wm}+“"'W-o 

where 


(38) 


a  gsdi*  1 
•*(r.  +  rfc) 

With  the  boundary  conditions  that  N(r)  is  symmetric  about  r  ■  0,  and 
smoothly  varying  so  that  ■  0  at  x  •  0,  the  solution  is 

*(')  -  Afo|l^*^ora]-*  (39) 

where  •  N(*)  at  x  •  0.  Equation  39  is  called  the  Bennett  distribution 
(Fig.  5). 


i 
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For  a  solid  state  material,  Na  can  be  related  to  the  intrinsic  and  in¬ 
jected  carrier  concentration.  Nj  will  have  a  temperature  dominated  factor  and 
an  injection  process  related  field  factor. 

To  relate  N(r)  back  to  Ni  of  Equation  25, 


X  f  2*r.S{r)<ir  m  J  2»r.V  ,(  1  -  .V  >,*)  ' 


•V  ( #)  -  V,rtf:(l  -  .V,6r:)  -■ 


But  the  plasma  extends  to  r  *  *  ,  so 


*.(*)  -  i 


If  a  is  defined  as  a  fraction  of  Hi  (oc)  in  a  finite  cylinder  of  radius. 


.¥,(/*)  -  «*,(<»)  »  o- 


substltutlng  what  Into  Equation  40  results  In 


(*)' 


If  a  *  0.9,  or  90  percent  of  the  carriers  are  within  the  cylinder  of 
concern,  the  last  equation  becomes 

sfN^kR  -  3 

defining  the  relation  of  R  to  our  problem  In  a  simple  manner. 
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The  Bennett  distribution  can  now  be  used  to  derive  the  applicable  P(r) 
and  B( r)  values 

J,(r)  *  .V(r)«|*.  ♦  V*  *  JV„|  1  ♦  -*■  W*1 

*  fn 

I ,(r)  »  J  2nrS,Al  ♦  |v.  ♦  *  2».v„e|v.  ♦  *h\  J  r(l  ♦  hS^r7)-3* r 

/,(-)  *  2w.S\,*{v.  *  V*)  (  —  (7,  k.\„R'-)  *  26,V,,( i  -  bS„r2)  ) 

'  . . .  *  ((I  -  »J»„rJ|(l*  »««»  =  )) 


for 


B.(r)  - 


i*o*\y.  +  vJ  (  9  -  No**3 


(  9  -  N0*r*  \ 
[{r  +  Notr*)) 


and  for  P(r) 


Substitution  of  fl,(r)  above  into  P(r)  Is  extremely  difficult  to  solve  and 
yields  the  following  relation  after  the  use  of  Integral  tables: 

9  /  1  1  R3  \ 

[T)  "  l2  {R^  +  NoiR3)  *  R^  +  NohR3)3  +  lo#(l  +  S^R3)  ) 

9/  1  1  r2  \ 

"  2  \r(l  ♦  Notr*)  *  r(l  +  AM'2)3  *  <0,(1  +  AM'3)  j 

/  AM  >\ , 

+  V  4*u4(  1  *  y<yiR3)3  '  4AM<  1  «■  AMr2)2  / 1 


where 


-AM3Moj*«  v*»la 


The  form  of  N(r)  In  this  example  possesses  attributes  of  the  strongly  and 
weakly  ionized  plasma.  Near  the  center  of  the  column,  the  plasma  is  strongly 


1* 
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ionized  and  slackens  in  activity  as  r  increases  in  distance  from  the 
center. 

6.3  THE  DYNAMIC  PINCH 

If  the  equilibrium  condition  for  the  equilibrium  pinch  cannot  be  met, 
the  dynamic  pinch  must  be  applied. 

By  starting  with  a  known  input  current 


and 


fl.(r) 


Mof(r) 

2wr 


Bj(r) 
2  MO 


N{r,t)kWr,t)  +  Tk{r,t)) 


can  be  used  to  find  the  inward  force  on  a  collapsing  current  sheath.  As  the 
sheath  collapses  inward,  it  collects  plasma  in  front  of  it.  Hence,  this  model 
is  called  the  snowplow  model.  The  collapsing  force  Is  equated  to  the  plasma's 
rate  of  change  of  momentum  to  arrive  at  the  equation  of  motion. 

In  a  manner  similar  to  a  mechanical  system,  the  snowplow  model  shows 
(Ref.  11)  that  a  shock  wave  Initiated  by  the  pinch  reflects  off  the  axis  and 
back  to  the  pinch  front  resulting  In  an  oscillation. 

Instead  of  displaying  a  constant  pinch  radius,  a  dynamic  pinch  plasma 
oscillates  around  the  equilibrium  radius  established  by  the  input  current  and 
plasma  parameters.  This  effect  has  been  called  bounce.  If  the  input  current 
oscillates,  it  is  possible  to  initiate  an  instability  that  will  continue 
throughout  current  application  (see  experimental  evidence  in  Ref.  10). 

Dynamic  pinch  problems  are  difficult  to  verify  when  applied  to  the  solid 
state  plasma.  Use  of  Monte  Carlo  modeling  should  begin  to  eliminate  problems 
in  modeling  the  dynamic  pinch. 
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7.0  GROWTH  INSTABILITIES  OF  THE  EQUILIBRIUM  SOLID  STATE  PLASMA 


As  mentioned  in  the  section  on  equilibrium  and  instabilities,  four 
different  categories  of  instabilities  are  possible:  Streaming  instabilities, 
Rayleigh-Taylor  instabilities.  Kinematic  instabilities  and  Universal  insta¬ 
bilities.  These  four  categories  can  be  placed  under  the  general  categories 
of  configuration  space  instabilities  and  velocity  space  instabilities. 

These  last  two  categories  are  useful  in  explaining  the  type  of  instabili 
ties  formed  in  semiconductor  plasma  formation  from  electrical  or  optical 
excitation. 


Under  electrical  excitation  the  transition  from  single-  to  double-carrier 
injection  under  low  fields  is  the  main  impetus  to  form  a  weakly  ionized  plasma. 
Avalanche  breakdown,  high  electric  fields,  or  high  temperatures  usually  form 
strongly  ionized  plasmas.  The  differences  between  these  two  conditions  also 
determine  the  type  of  instability  that  can  occur. 

7.1  SAUSAGE  INSTABILITY 


¥ 


If  an  otherwise  cylindrical  plasma  is  perturbed  by  an  external  force 
causing  a  spatial  variation  of  the  plasma,  a  new  equilibrium  position  will  be 
reached.  But,  in  moving  to  a  new  equilibrium  position,  free  energy  may  be 
released  resulting  in  an  instability. 

As  Fig.  6  shows,  bulges  can  form  in  the  plasma  as  stationary  or  moving 
waves.  The  increased  current  density  in  pinched  portions  of  the  plasma  will 
create  an  Increased  magnetic  field  around  the  pinch.  In  areas  where  the  plasma 
has  extended,  the  decreased  current  density  lowers  the  magnetic  field  strength 
allowing  plasma  pressure  to  expand  outward.  These  mechanisms  can  work  in 
concert  to  actually  pinch  the  column  off.  Externally,  rapid  oscillations  in 
current  are  measured  as  this  instability  occurs. 

In  a  semiconductor  diode  (PN  or  PIN),  this  type  of  instability  readily 
occurs  due  to  the  depletion  layer  configuration.  Since  current  streams  out  of 
the  plasma  column  ends,  and  kinetic  and  magnetic  pressures  are  highest  at  the 
center  of  the  column  length,  the  sausage  Instability  forms  upon  the  plasma 
reaching  the  depletion  edge  (Fig.  7). 

Similar  situations  can  occur  for  PIN  diodes,  and  transistors  for  collector- 
emitter,  collector-base  and  emitter-base  plasma  formation.  Avalanche  breakdown 
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the  instability  grows,  B*  will  almost  pinch  off  current  at  troughs,  ani 
kinetic  pressure  will  elongate  crests  at  B#  .  Current  measured  during 
this  instability  will  be  sinusoidal  in  nature. 


AFWL-TR-85-116 
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Figure  7.  Depiction  of  possible  areas  of  plasma  formation  in 
PN  and  PIN  diodes.  The  shaded  areas  show  a  plasma 
column  with  fluted  ends.  The  shape  is  indicative  of 
a  sausage  instability  which  may  occur  upon  pinch  due 
to  the  devices  geometry.  Most  experiments  that  study 
plasma  pinch  in  semiconductors  see  oscillations  simul¬ 
taneously  with  plasma  formation. 
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in  the  midst  of  plasma  formation  will  cause  a  perturbation  in  the  column  that 
will  resemble  the  sausage  instability. 

If  avalanche  has  already  occurred,  initiating  plasma  formation,  extensive 
heating  near  the  pinch  center  could  cause  a  sausage  instability.  This  event 
could  also  be  explained  as  the  beginning  of  a  thermal  pinch,  to  be  descried 
later  in  the  report. 

7.2  KINK  INSTABILITY 

If  the  plasma  column  is  disturbed  by  a  relatively  large  perturbation 
caused  by  device  construction,  the  kink  instability  occurs. 

As  Fig.  8  illustrates,  the  magnetic  field  lines  are  more  concentrated  in 
Region  1  than  in  Region  2.  Thus,  the  magnetic  pressure  is  greater  in  Region  1 
than  in  Region  2.  A  kink  started  in  this  manner  uses  this  feedback  mechanism 
to  easily  continue  the  initial  perturbation. 

Because  of  transistor  construction  using  mesa  architecture  (as  shown  in 
Fig.  9),  a  kink  instability  due  to  base  current  flow  can  easily  form. 

If  a  surface  is  near  the  direction  of  Instability,  surface  recombination 
could  aid  instability  growth.  In  a  transistor,  such  action  would  appear  to 
increase  the  base  resistance  to  current  flow  and  debias  the  emitter  under 
collector-emitter  breakdown  in  the  common-emitter  mode. 

A  second  cause  for  the  kink  instability  is  plasma  formation  in  depletion 
regions  that  curved  perpendicularly  to  current  flow. 

7.3  STREAMING  INSTABILITIES 

It  is  well  known  from  classical  mechanics  that  two  dynamic  systems  inter¬ 
act  most  effectively  if  the  system  velocities  or  frequency  of  oscillations 
are  about  the  same. 

In  a  semiconductor,  the  mass  difference  between  holes  and  electrons  pro¬ 
duces  two  separate  dynamic  systems  for  the  semiconductor  plasma.  The  plasma, 
as  a  whole,  oscillates  at  a  frequency  (separate  from  the  electron  and  hole 
plasma  frequencies)  in  such  a  way  that  its  phase  velocity  is  normally  higher 
than  the  carriers'  thermal  velocity  by  a  large  magnitude,  thereby  preventing 
any  exchange  of  energy. 
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Figure  8.  Kink  instability.  The  perturbed  column  creates  a 
curvature  that  increases  the  magnetic  pressure  in 
region  (1)  at  the  expense  of  lessening  it  in  region 
(2).  The  kinetic  pressure  responds  by  exerting 
greater  pressure  towards  region  (2).  The  feedback 
forces  the  kink  to  continue  to  grow  in  the  direc¬ 
tion  of  region  (2). 
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Figure  9.  A  mesa  P-N-P  transistor  in  the  common-emitter  mode 
in  cut-off.  As  current  is  withdrawn  from  the  base, 
plasmas  may  form  and  current  crowding  may  instigate 
kink  instabilities  leading  to  further  current  crowding. 
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As  carriers  heat  up,  their  velocities  approach  the  plasma's  phase 
velocity.  As  their  velocities  surpass  the  phase  velocity,  the  plasma  oscilla¬ 
tions  gain  energy  from  their  carriers'  velocities,  slowing  some  carriers  down 
and  speeding  others  up.  This  process  changes  the  normal  one-hump  Maxwellian 
velocity  curve  to  a  two-hump  curve.  The  plasma  oscillations  grow  in  magnitude 
from  the  external  energy  source  resulting  in  the  two  stream  instability 
(Ref.  12,  Fig.  10). 

The  growing  plasma  oscillations  begin  to  transfer  energy  back  to  the 
lattice  via  optical  or  acoustic  phonon  excitation,  heating  the  lattice  till 
it  melts  and  enters  the  thermal  pinch  phase. 


The  two-stream  instability  can  occur  via  excitation  of  the  plasma  by 
heating  or  an  increasing  electric  field.  It  can  also  occur  as  a  consequence 
of  the  quantum  mechanical  aspects  of  a  semiconductor  material. 


The  formula  relating  the  carrier  velocity  to  plasma  phase  velocity  for 
occurrence  of  the  two-stream  instability  was  discerned  by  Pines  and  Schrieffer 
(Ref.  13)  as  the  following  in  the  long  wavelength  limit. 


If  r€-rfc  ,  then  the  threshold  is  /  /  \  i\ 

—«■*(*♦  (5)) 

when  the  relative  drift  velocity,  V  ,  is  equal  to  the  electron  thermal 
velocity. 


However,  mg  and  m^  are  dependent  on  cyrstal  orientation  and  carrier 
excitation  level.  As  carriers  gain  enough  energy  to  enter  the  conduction 
band,  they  eventually  enter  other  minima  that  reduce  the  carrier's  effective 
mass;  thus,  increasing  the  value  of  Vt  and  «D,  which  interact  with  higher 
energy  phonons  in  a  strongly  coherent  manner. 


This  type  of  plasma  formation  would  occur  for  very  rapid  excitation  of  a 
semiconductor  by  an  optical  or  electric  field,  and  should  be  dependent  on  the 
crystal  orientation  and  band  structure.  The  occurrence  of  this  instability 
would  rapidly  degrade  to  a  thermal  pinch  under  room  temperature  conditions. 
Actual  measurement  of  this  instability  has  been  accomplished  at  low  tempera¬ 
tures  for  n-type  InSb  (Ref.  14). 
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Figure  10.  Unstable  Maxwellian  distribution.  The 
unstable  distribution  is  produced  by  the 
two-stream  instability. 
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8.0  THERMAL  PINCH 

At  the  high  temperatures  produced  by  semiconductor  plasmas  under  extreme 
conditions,  a  type  of  instability  occurs  that  is  the  result  of  a  magneto-thermal 
pinch  transitioning  to  a  pinch  dominated  by  thermodynamic  driving  forces. 

This  type  of  pinch  has  been  labeled  the  thermal  pinch  (Ref.  15)  to  dis¬ 
tinguish  it  from  the  previous  discussion  of  a  magneto-thermal  pinch. 

When  heating  in  the  core  of  a  magneto- thermal  pinch  reaches  levels  where 
the  intrinsic  carrier  concentration  rises  rapidly  to  create  a  degenerate 
plasma,  the  kinetic  pressure  quickly  overcomes  the  magnetic  pressure.  The 
material  begins  to  change  phase  from  a  solid  to  a  liquid  and  expands  rapidly 
to  engulf  the  available  region.  This  type  of  pinch  can  literally  explode  out 
the  ends  of  the  semiconductor  being  pulsed  and  totally  dominates  the  prefailure 
electrical  characteristics. 

Proof  for  this  type  of  pinch  was  discovered  by  Betsy  Ancker-Johnson  and 
James  Drummond  In  References  10  and  15.  Their  experiments  with  InSb  actually 
formed  hollow  melt  channels  In  the  specimens  examined.  Additional  measurements 
showed  a  distinct  trend  for  pinch  channel  radius  found  via  Bennett's  tempera¬ 
ture  (using  the  magneto- thermal  pinch)  to  move  toward  a  channel  radius  dictated 
by  the  semiconductor  melting  temperature  (or  thermal  pinch)  as  shown  in  Fig.  11. 

Since  the  thermal  pinch  literally  moves  through  a  material  In  a  manner 
similar  to  electromigration,  its  velocity  of  propagation  is  highly  dependent 
on  the  material  configuration  and  thermodynamic  properties  (heat  of  fusion, 
coefficient  of  expansion,  available  energy  to  do  work  in  the  plasma,  and 
geometry) . 

The  dependence  of  the  thermal  pinch  on  so  many  micro-  and  macroscopic 
variables  makes  modeling  of  a  thermal  pinch  in  devices  very  difficult,  and 
somewhat  probabilistic. 

This  is  the  type  of  plasma  formation  that  occurs  in  semiconductor  devices 
pulsed  to  destruction.  The  character  of  the  thermal  plasma  pinch  would  add  a 
variable  character  that  may  explain  the  variation  found  in  experimental  studies 
(Ref.  1). 
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Figure  11.  The  pinch  radius  divided  by  the  plasma  current 
as  a  function  of  power  input  per  unit  length  of 
plasma:  Experimental  data  at  the  onset  of 
pinching  (x),  at  the  end  of  the  pulses  (o).  The 
upper  curve  is  the  result  of  magnetic  theory  and 
the  lower,  thermal  pinch  theory.  The  numbers 
associated  with  the  arrows  correspond  to  the  net 
energy  transferred  to  the  pinch  column  in  units 
of  energy  required  to  melt  the  column  (Ref.  15). 
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Theoretical  analysis  of  device  failure  based  on  thermal  models  would  then 
be  able  to  .model  some  general  trends,  but  would  be  totally  inadequate  to  ex¬ 
plain  propagation  of  the  thermal  pinch  or  the  initiating  conditions  prompted 
by  the  magneto- thermal  pinch.  These  models  do  not  simulate  plasma  effects  or 
changing  solid  state  parameters  that  affect  the  plasma. 

In  optical  induced  failure,  the  thermal  pinch  also  represents  an  extreme 
of  the  pinch  phenomena.  Once  thermal  pinching  occurs,  the  plasma  would 
blossom  outward  actually  melting  the  material.  The  model  put  forth  by  Ancker- 
Johnson  et  al  in  References  15  and  16  should  provide  a  good  starting  point  to 
determine  the  initial  conditions  surrounding  optically  induced  failure  in 
semiconductors. 

The  Monte  Carlo  program  produced  for  this  effort  still  must  be  modified 
to  model  the  thermal  pinch. 
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9.0  MONTE  CARLO  METHOD 


Damage  to  semiconductor  devices  by  intense  electromagnetic  pulses  has  been 
studied  over  many  decades.  However,  understanding  of  damage  mechanisms  at 
the  level  of  basic  physical  processes  is  still  very  incomplete.  This  is  not 
surprising  in  view  of  the  complexity,  in  terms  of  solid-state  device  physics, 
of  the  phenomena  involved.  With  the  intense  fields  that  cause  damage,  basic 
analysis  requires  much  more  accurate  energy  band  structure  and  carrier  scat¬ 
tering  rates  than  have  been  adequate  heretofore  even  for  most  problems  in  high 
field  transport.  Moreover,  damage  occurs  typically  through  a  series  of  rather 
distinct  stages.  A  complete  theoretical  treatment  also  involves  modeling 
collective  effects,  such  as  space  charge  and  carrier-carrier  scattering,  in 
order  to  describe  the  later  stages.  However,  it  can  now  be  said  that  all  of 
these  difficulties  are  tractable.  Spurred  by  Its  importance  in  microelectron¬ 
ics,  progress  in  the  numerical  calculation  of  high-field  transport  effects  has 
produced  computational  techniques  for  treating  these  necessary  refinements 
(Refs.  17  and  18).  It  remains  now  to  combine  the  various  methods  and  apply 
the  resulting  computational  power  to  the  problem  of  device  damage. 

A  damage  mechanism,  such  as  second  breakdown  (Ref.  19),  occurs  in  a 
device  when  the  electromagnetic  field  exceeds  a  threshold  value.  It  evolves 
from  an  Initial  precipitative  stage  of  carrier  growth,  through  one  or  more 
stages  where  the  effects  of  large  current  and  charge  densities  are  important, 
to  a  final  phase  of  material  breakdown.  In  general ,  space,  as  well  as  time 
dependence,  enters  (particularly  in  the  later  stages  of  plasma  development) 
even  for  the  case  of  damage  in  the  homogeneous  bulk. 

Full  understanding  of  a  phenomenon  requires  its  description  at  the  micro¬ 
scopic  level.  Consider  the  mechanism  by  which  negative  differential  mobility 
occurs  in  the  bulk  of  certain  semiconductors.  This  hot  electron  phenomenon 
is  a  transferred  electron  effect  (Ref.  19).  When  the  electric  field  exceeds 
some  threshold  value,  electron  scattering  from  the  lowest  valley  to  upper 
valleys  of  the  conduction  band  becomes  large  enough  to  cause  significant 
repopulation  to  the  upper  valleys.  The  (steady  state)  mobility  thereafter 
decreases  with  increasing  field.  Once  this  basic  description  is  known,  one 
has  a  unifying  theory  whereby  it  is  possible  to  predict  the  occurrence  and 
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magni  tude  of  the  effect  in  various  semiconductors  from  a  knowledge  of  their 
band  structure  and  scattering  properties. 

When  an  external  electromagnetic  field  is  applied  to  a  conducting 
material,  the  most  complete  quasiclassical  description  of  the  resulting  charge 
transport  is  given  by  the  carrier  distribution  function  f(r,  p,  t).  This 
specifies  the  statistical  distribution  of  mobile  charge  carriers  over  phase 
(momentum  p  plus  coordinate  r)  space.  In  the  absence  of  both  external  fields 
and  temperature  gradients,  the  distribution  function  is  everywhere  independent 
of  time  and  is  that  of  thermal  equilibrium,  viz.,  a  Fermi-Dirac  distribution. 
Upon  application  of  a  field  which  is  thereafter  time  independent  dc  or  perio¬ 
dic  of  constant  amplitude,  the  distribution  rapidly  attains  a  steady  state 
form  (possibly  with  spatial  variation)  that  describes  the  steady-state  charge 
transport  resulting.  In  a  truly  transient  situation,  such  as  that  of  a  damage 
mechanism,  the  distribution  function  will  have  temporal  as  well  as  possible 
spatial  dependence  as  the  carriers  respond  change  in  the  field  acting  on  them. 


In  the  quasi-classlcal  picture  of  charge  transport,  the  mobile  carriers 
are  approximately  Independent  particles.  For  crystalline  materials,  carriers 
occupy  the  stationary  quantum  Bloch  states  of  the  perfect  ionic  lattice.  The 
free-electron-1 ike  Bloch  states  are  each  specified  by  a  value  of  crystal- 
momentum  vector  p,  and  the  energy  eigenvalue  of  a  state  is  given  by  the 
applicable  band  structure  function  E(]>).  The  state  of  a  carrier  will  change 
as  the  result  of  deviations  from  perfect  crystal  periodicity  caused  by  phonons 
(lattice  vibrations),  impurities  and  imperfections,  and  by  the  presence  of 
other  carriers.  The  transition  of  a  carrier  state  that  occurs  in  an  inter¬ 
action  with  one  of  these  perturbation  sources  is  described  as  a  scattering 
(or  collision)  event.  A  scattering  event  is  governed  by  the  laws  of  quantum 
mechanics.  For  each  type  of  scattering  process,  there  is  a  probabilistic 
scattering  function  ("p,  c;  p',  c"),  which  is  the  transition  probability  per 
unit  time  for  a  carrier  in  state  "p  with  the  crystal  in  state  c  to  scatter  into 
the  state  ~p'  leaving  the  crystal  in  state  c". 

In  the  thermal  equilibrium,  the  carrier  distribution  gives  detailed  bal¬ 
ance  between  any  scattering  event  and  its  time-reversed  event.  The  thermal 
motion  of  the  carriers  is  completely  chaotic  with  no  net  transport  of  charge 
or  energy.  Where  no  built-in  field  is  present,  the  Bloch  state  of  a  carrier 
remains  unchanged  between  collisions. 
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Charge  transport  occurs  if  an  external  field  is  applied.  In  free  flight 
between  two  consecutive  scatterings,  a  carrier  moves  solely  under  the  force 
of  the  electromagnet ic  field  acting  on  it  and  its  momentum  changes  in  accor¬ 
dance  with  the  classical  laws  of  motion.  Each  free  flight  is  randomly  inter¬ 
rupted  when  the  carrier  interacts  with  the  lattice  or  with  another  carrier. 
Thus,  the  state  of  a  carrier  varies  continuously  through  the  Bloch  states  as 
it  is  accelerated  by  the  field  during  a  free  flight.  But  the  state  changes 
abruptly  and  the  change  can  be  large  in  the  scattering  event  that  terminates 
the  flight.  The  distribution  function  in  this  quasi -classical  picture  is 
governed  by  the  Boltzman  transport  equation  (BTE),  which  incorporates  the 
effects  on  the  distribution  of  the  dynamical  processes  (carrier  acceleration 
by  the  field  and  carrier  collisions)  and  of  spatial  or  temporal  variations 
in  the  distribution.  The  departure  of  the  distribution  function  from  equili¬ 
brium  with  the  lattice  increases  with  the  strength  of  the  applied  field.  It 
is  usually  permissible  to  consider  the  phonon  distribution  as  always 
thermal i zed . 

The  values  of  the  distribution  function  at  relatively  fine-grained  cells 
in  phase  space  comprise  a  huge  amount  of  information,  often  more  than  is 
desired.  Physical  quantities  characterizing  carrier  motion,  such  as  carrier 
energy  E(p)  and  velocity  V^E(p),  are  functions  of  momentum.  Significant 
parameters  in  charge  transport  are  the  average  of  such  quantities  over  the 
distribution  function  throughout  momentum  space  or  in  a  certain  region  of 
momentum  space  (such  as  an  energy  valley).  In  most  cases,  the  salient  features 
of  a  phenomenon  are  revealed  by  the  behavior  of  a  relatively  few  transport 
parameters  and  by  the  coarse-grained  dependence  of  the  distribution  function 
on  some  momentum  component  or  on  energy. 

Detailed  analysis  of  damage  at  the  level  of  fundamental  physical  processes 
'■equires  more  powerful  and  complex  computational  methods  than  have  previously 
been  brought  to  bear  on  the  problem.  But  numerical  analysis  of  high-field 
(hot-carrier)  behavior  of  semiconductor  devices  has  now  advanced  to  where 
such  first-principles  modeling  can  be  applied  to  the  extreme  fields  that 
cause  damage.  Once  the  roles  of  basic  physical  properties  are  properly 
understood  for  a  particular  mechanism,  a  reliable  theoretical  basis  will  be 
available  for  developing  devices  more  resistant  to  damage  from  this  cause. 

One  can  hope  thereby  to  make  optimum  decisions  on  material  selection,  device 
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design,  and  circuit  characteristics.  An  additional  benefit  will  accrue  in 
the  area  of  device  engineering.  The  simplified  equations  commonly  used  in 
device  analysis  and  design  can  be  made  more  accurate  and  confidence  limits 
set  for  their  application  in  the  damage  regime.  The  transport  parameters 
appearing  in  these  equations  can  be  calculated  with  the  help  of  detailed 
numerical  methods  and  their  predictions  can  be  compared  with  the  relatively 
exact  calculations  (Ref.  20). 

9.1  NUMERICAL  ANALYSIS  OF  SEMICONDUCTOR  DEVICES 

To  properly  understand  and  design  today's  electronic  devices  requires 
detailed  modeling  of  the  space-  and  time-dependent  charge  and  current  distri¬ 
butions  at  the  operating  potentials  of  interest.  Numerical  computations  are 
generally  resorted  to  (1)  when  the  device  equations  do  not  admit  analytic 
solutions  because  of  complicated  transient  effects  or  spatial  structures, 

(2)  when  device  transport  is  nonlinear  because  of  hot-carrier  contributions, 
and  (3)  in  computer-aided  design  of  devices  that  require  multidimensional 
structures  for  their  operation. 

The  behavior  of  a  solid  state  device  whose  dimensions  are  greater  than 
one  micron  has  conventionally  been  analyzed  (particularly  in  the  case  of 
silicon)  by  means  of  the  semiconductor  equations.  These  equations  are  based 
on  the  classical  nydrodynamical  model  of  carrier  current,  the  equation  for 
the  sel f-consistent  potential  in  the  device.  Current  continuity  is  obtained 
as  the  zeroth  momentum  moment  of  the  BTE,  while  the  classical  expression  for 
carrier  current  represents  the  drift-diffusion  approximation  to  the  momentum 
balance  equation.  Momentum  balance  is  obtained  as  the  first  momentum  moment 
of  the  BTE.  Poisson's  equation  must  be  solved  self-consistently  with  the 
prescribed  potentials  applied  at  the  boundaries  and  the  charge  density  given 
by  the  other  two  equations. 

In  the  high-field  regime  where  the  current  response  is  nonlinear,  the 
mobility  and  diffusion  coefficient  are  both  field  dependent  and  must  be  known 
from  laboratory  measurements  or  from  exact  calculations.  Recent  studies  have 
compared  high-field  transport  behavior  predicted  by  the  semiconductor  equa¬ 
tions  with  more  accurate  solutions  of  the  BTE  (Ref.  19).  These  show  that 
the  conventional  equations  have  limited  utility  in  submicron  and  in  high- 
frequency  devices  and  are  generally  inaccurate  in  high-field  transient 
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situations.  In  the  classical  current  model,  it  is  assumed  that  (1)  carrier 
current  responds  instantaneously  to  changes  in  electric  field,  (2)  the  cor¬ 
responding  steady-state  mobility  and  diffusion  coefficient  can  be  used  at 
each  field  value,  and  (4)  diffusion  and  drift  do  not  influence  each  other. 
However,  fast  transient  response  is  nonlocal  in  time  and  requires  the  BTE 
for  an  accurate  description.  The  semiconductor  equations  are  inherently 
unable  to  describe  transient  phenomena  occurring  on  a  time  scale  of  pico¬ 
seconds  or  less,  as  for  example  velocity  overshoot. 

A  set  of  exact  moment  (or  balance)  equations  from  the  BTE  must  be  used 
first  to  obtain  more  accurate  device  equations.  In  one  approach  the  electron 
distribution  function  for  each  energy  valley  is  approximated  by  an  analytic 
form,  e.g.,  a  Maxwellian  distribution  displaced  in  momentum  space  so  that 
the  average  carrier  momentum  pa  is  along  the  electric  field  direction.  The 
analytic  form  depends  on  a  small  number  of  time-dependent  parameters  (e.g., 

Pd  electron  temperature,  and  concentration)  which  are  substituted  into  a 
like  number  of  balance  equations  for  the  particular  valley  (starting  with  the 
zeroth  moment  equation).  The  resulting  set  of  coupled  equations,  containing 
the  parameters  as  unknowns,  is  used  to  model  device  behavior.  Scattering 
integrals  appearing  in  the  balance  equations  may  be  given  in  approximate 
analytic  form  or  evaluated  numerically. 


When  the  drifted-Maxwel 1 ian  distribution  is  not  justified,  the  semi- 
empirical  transport  equation  method  (SETEM)  can  be  employed.  In  this  approach 
one  starts  with  the  balance  equations  for  charge,  momentum  and  energy  along 
with  Poisson's  equation.  The  scattering  integrals  in  the  momentum  and  energy 
balance  equations  are  approximated  with  effective  momentum  and  energy  relax¬ 
ation  rates,  respectively,  which  are  field  dependent.  At  each  field  value 
of  interest,  the  relaxation  rates  are  chosen  to  satisfy  the  steady-state 
form  of  the  balance  equations  with  exact  values  of  steady-state  drift  velo¬ 
city,  average  energy  and  average  effective  mass  inserted  in  the  equations. 

The  exact  values  required  are  obtained  from  detailed  numerical  solutions  of 
the  steady-state  BTEs,  usually  by  Monte  Carlo  simulation.  Thus,  the  SETEM 
attempts  to  describe  nonstationary  behavior  with  the  steady-state  parameters. 

The  accuracy  of  the  more  elaborate  device  equations  just  described 
ranges  from  very  good  to  exceedingly  poor,  depending  on  the  validity  of  the 
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approximations  used  in  the  particular  situation  treated.  Surprisingly  good 
results  can  be  obtained,  but  these  approaches  are  applicable  in  a  limited 
number  of  cases. 

The  most  accurate  description  of  high-field  transport  behavior  that  is 
computationally  tractable  is  the  Boltzmann  transport  equation.  Accurate 
numerical  calculations  with  the  BTE  retain  in  effect  the  moment  equations  to 
all  orders  and  can  describe  quite  general  temporal  and  spatial  effects. 

Various  numerical  techniques  for  the  solution  of  the  BTE  have  been  developed. 
Their  use  has  been  limited  in  the  past  mainly  by  the  computational  expense 
involved.  But  with  the  rapid  development  of  fast  large-memory  computers,  the 
cost  factor  is  much  less  important  now  and  should  soon  be  a  relatively  minor 
consideration. 

Of  the  various  numerical  approaches  to  the  BTE,  the  Monte  Carlo  method 
is  now  by  far  the  most  popular  in  solid  state  transport  analysis.  It  is  easy 
to  program  and  to  include  realistic  band  structure  and  scattering  models. 
Extension  to  time-  and  space-dependent  phenomena  is  readily  implemented  by  an 
ensemble  simulation.  It  is  the  only  method  by  which  fluctuation  phenomena 
can  be  studied  with  a  straightforward  analysis.  Being  a-  direct  simulative 
procedure,  it  is  the  most  directly  interpretable  from  the  physical  point  of 
view.  The  analysis  of  any  phenomenon  on  the  microscopic  level,  i.e.,  in  terms 
of  band  structure  and  scattering  processes,  is  quite  transparent.  Simulation 
with  the  self-consistent  field  can  be  performed  in  a  direct  way  in  ensemble 
Monte  Carlo.  Inclusion  of  short-range  carrier-carrier  interactions,  however, 
requires  an  elaboration  of  Monte  Carlo  that  is  more  difficult  to  implement 
but  straightforward  in  principle  (Ref.  22).  Because  of  its  comparative  advan¬ 
tage  over  other  detailed  numerical  techniques,  the  application  of  Monte  Carlo 
simulation  to  device  analysis  and  design  is  increasing  steadily  (Ref.  17). 

9.2  QUASI-CLASSICAL  DESCRIPTION  OF  CHARGE  TRANSPORT 

When  treating  the  interaction  of  an  applied  field  with  a  material,  the 
basic  formulation  of  charge  transport  is  almost  always  quasi-classical  in 
nature.  The  mobile  charge  carriers  in  semiconductors  are  electrons  and  holes. 
The  electron  and  hole  systems  are  given  a  quasi-classical  description  in  which 
an  individual  mobile  carrier  is  assigned  a  Bloch  energy  band  n  and  a  wave- 
vector  k,  specifying  the  Bloch  energy  state,  and  is  simultaneously  assigned 


a  position  r.  The  Bloch  state  has  associated  crystal  momentum,  p  =  hk  ,  and 
energy,  EnCk). 

The  81och  states  are  those  of  a  perfect  crystal  lattice.  As  a  carrier 
moves  through  the  lattice,  its  Bloch  state  is  randomly  perturbed  by  crystal 
imperfections,  namely,  lattice  vibrations  (or  phonons),  impurities,  structural 
defects,  interfaces,  and  by  other  carriers.  The  sudden  state  changes  are 
pictured  as  various  scattering  processes  in  which  transition  to  another  Bloch 
state  occurs  by  carrier  interaction  with  one  of  the  perturbing  sources.  The 
scattering  processes  obey  the  laws  of  quantum  mechanics  and  their  cross  sec¬ 
tions  are  derived  with  the  quantum  theory  of  scattering. 

In  the  quasi-classical  description,  the  state  of  the  carrier  system  at 
time  t  is  specified  by  a  distribution  function  f p ("k »  7,  t)  for  each  relevant 
energy  band  n.  The  number  of  carriers  in  the  band  at  time  t  occupying  the 
volume  elements  d37  at  ¥  and  d3r  at  r  is  (2it)"3fn  ("k,  7,  t)  d3  k  d1  r  where  n 
denotes  the  band  (or  the  particular  minimum  of  a  band  and  possibly  the  spin 

A 

component  value  as  well)  and  f  is  the  absolute  distribution  function.  During 
free  flight,  i.e.,  on  a  trajectory  connecting  two  successive  scattering  events, 
the  carrier  motion  is  governed  by  the  field  acting  on  it  and  obeys  the  quasi- 
classical  dynamical  laws: 

(*4) 

k.P  (45) 

Here,  vjr  is  the  carrier  velocity  in  the  state  7  and  7  is  the  (self- 
consistent)  force  field  (apart  from  the  periodic  field  of  the  perfect  ionic 
lattice)  in  which  the  carrier  moves  between  collisions. 

When  there  is  no  applied  field  and  no  temperature  variation,  the  carriers 
are  in  equilibrium  with  the  lattice  and  fn(7,  7)  is  everywhere  a  Fermi -Dirac 
distribution  with  temperature  equal  to  that  of  the  lattice.  In  response  to 
an  external  field  the  distribution  changes.  From  the  nonequilibrium  distri¬ 
bution  function  at  time  t,  all  transport  parameters  (mean  carrier  energy, 
drift  velocity,  average  collision  rates,  etc.)  can  be  calculated  at  t. 

The  most  general  transport  theory  within  the  context  of  the  quasi- 
classical  description  is  essentially  that  embodied  in  the  set  of  coupled 
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Boltzmann  equations  for  the  distribution  functions.  Of  the  various  computa¬ 
tional  methods  for  analyzing  charge  transport  in  semiconductors  with  the 
full  accuracy  of  the  Boltzmann  equation,  Monte  Carlo  simulation  has  emerged 
as  the  most  useful  technique  in  high-field  problems.  Monte  Carlo  combines 
computational  simplicity  with  the  broadest  applicability  to  complex  band 
structures  and  scattering  processes,  and  is  also  the  most  powerful  method 
for  treating  problems  having  space  dependence. 

9.3  MONTE  CARLO  SIMULATION 

In  direct  Monte  Carlo  (Ref.  21)  simulation,  the  behavior  of  the  physical 
system  is  simulated  by  computing  individual  particle  motions.  The  calculated 
particle  histories  are  treated  as  a  sample  set  from  an  ensemble  of  histories 
for  the  particles  of  the  physical  system.  The  average  of  a  particle  variable 
over  this  set  at  a  conmon  time  is  an  estimator  for  the  average  over  the  phys¬ 
ical  ensemble  at  that  time. 

For  charge  transport  in  a  semiconductor,  the  particle  (electron  or  hole) 
history  (in  the  quasi-classical  approximation)  is  a  sequence  of  free  flights 
under  the  influence  of  the  electromagnetic  field  alternating  with  scatterings. 
Each  scattering  marks  the  termination  of  one  free  flight  and  the  start  of  the 
next  one. 

Clearly  a  Monte  Carlo  computer  program  for  charge  transport  calculations 
must  contain  algorithms  to  perform  the  following  basic  functions: 

1.  Describe  the  energy  band  structure,  En(k),  over  the  full  range  of 
particle  energies  that  are  important  in  the  problem. 

2.  Account  for  the  effect  of  the  applied  field  on  particle  motion  and 
energy  absorbed  from  the  field. 

3.  Calculate  the  scattering  function  and  scattering  rate  for  each  colli 
sion  process  as  required  in  the  simulation. 

4.  Statistically  select,  according  to  the  scattering  probability  func¬ 
tions  of  the  various  collision  processes,  values  of  the  random  quantities 
that  enter  into  the  motion  of  a  particle,  namely,  the  time  duration  of  a  free 
flight,  and  the  state  of  the  particle  immediately  after  the  scattering  event 
when  the  next  free  flight  commences. 
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A  particle  simulation  is  started  by  an  appropriate  statistical  selection 
of  the  initial  position  -m  phase  space  and  is  then  allowed  to  proceed  until 
the  particle  history  is  long  enough  (1)  to  include  a  number  of  scattering 
events  that  is  statistically  adequate  for  the  calculation,  or  (2)  to  traverse 
the  desired  space  or  time  interval  when  position  or  time  is  included  in  the 
particle  variables.  All  the  information  necessary  to  calculate  the  transport 
quantities  of  interest  (such  as  the  mean  particle  energy  or  mean  velocity) 
must  be  extracted  and  recorded  from  each  sinqle  free  fliqht. 

In  a  steady-state  situation  with  spatial  uniformity,  a  single  particle 
history  taken  over  a  sufficiently  long  time  can  provide  the  data  for  calculat¬ 
ing  the  transport  behavior.  For  such  steady-state  problems,  ergodicity  implies 
that  any  single  long  history  is  representative  of  all  the  particle  histories 
in  the  system.  Then  time-averaging  a  variable  over  the  particle  history  is 
equivalent  to  an  instantaneous  average  of  the  variable  over  the  distribution 
function  in  7-space. 

The  general  approach  of  using  time  averages  over  a  single  particle 
history  is  called  Single-Particle  Monte  Carlo  (SPMC).  The  time  spent  by  the 
particle  in  part  of  7-space  is  proportional  to  the  distribution  function. 
Setting  up  a  mesh  of  cells  in  7-space,  a  histogram  may  be  computed  which 
represents  the  time  spent  by  the  particle  inside  each  cell  of  the  mesh.  If 
a  mesh  with  sufficiently  fine  subdivisions  and  a  simulation  time  long  enough 
to  ensure  a  convergent  result  is  used,  the  histogram  will  accurately  estimate 
the  particle  distribution  function  and  can  be  used  to  calculate  mean  values. 

Usually  an  average  can  be  computed  more  efficiently  by  recording  the  mean 
value  for  each  free  flight,  accumulating  a  sufficient  number  of  such  flights, 
and  then  averaging  over  all  flights.  This  saves  computer  time  by  avoiding 
the  histogram  calculation  on  a  mesh  sufficiently  fine  to  ensure  accuracy  in 
averaging.  Other  estimators  for  average  values  in  steady  state  are  given  in 
the  literature  and  can  be  more  efficient  than  the  direct  cumulative  average 
over  flights  (Ref.  22). 

9.4  ENSEMBLE  MONTE  CARLO 

In  the  case  of  time  dependence  that  is  neither  periodic  nor  a  simple 
exponential  growth  or  decay  and  in  most  situations  of  spatial  dependence,  the 
simulated  history  of  a  single  carrier  does  not  adequately  represent  the 
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physical  ensemble  of  carriers.  The  general  approach  then  is  to  compute  a 
large  number  N  of  particle  histories,  which  are  treated  as  a  sample  set  from 
an  ensemble  of  particle  histories  that  corresponds  to  the  actual  physical 
system.  The  number  N  will  usually  be  much  smaller  than  the  number  of  carriers 
in  the  physical  system  being  analyzed.  It  need  only  be  large  enough  for  the 
system  quantities  of  interest  to  be  calculated  accurately  and  without  undue 
fluctuations  as  N  is  varied.  The  required  value  of  N  will  depend  on  the 
particular  problem. 

To  allow  for  interaction  between  articles  (see  Collective  Effects,  Section 
9.10),  the  particle  histories  are  computed  in  parallel.  For  each  particle, 
the  initial  Bloch  state  (and  position  also  when  there  is  spatial  nonuniformity 
initially)  is  selected  by  sampling  the  initial  distribution  of  the  physical 
system.  For  each  particle,  the  values  of  the  set  of  recorded  particle  vari¬ 
ables  at  a  designated  observation  time  are  used  to  estimate  the  desired 
transport  quantities  from  appropriate  averages  over  the  N  particles.  This 
gives  estimates  of  these  transport  parameters  (with  possible  spatial  variation) 
at  the  average  time  associated  with  a  set  of  observation  times.  In  the  next 
cycle  of  the  computation,  an  increment  of  each  particle  history  is  generated 
corresponding  to  a  time  step. 

The  evolution  of  this  simulation  at  successive  time  steps  represents  the 
actual  time  dependence.  This  is  referred  to  as  the  Many-Particle  Ensemble 
Monte  Carlo  (MPEMC)  scheme.  For  combined  time  and  space  dependence,  i.e., 
transient  dynamic  response  that  is  dissimilar  at  different  space  points, 

MPEMC  offers  the  only  tractable  approach  to  the  Boltzmann  equation. 

The  true  time  scheme  has  been  adopted  as  the  procedure  for  implementing 
the  program  cycles  that  generate  the  time  steps.  In  this  procedure,  the 

particle  histories  are  computed  between  successive  observation  times,  t^, 

(2) 

V  ,  ...  that  are  predetermined  and  are  the  same  for  each  particle.  The  time 
steps  At  =  t  +  ^  -  t^  can  all  be  equal,  or  can  be  spaced  more  closely 


✓ , 


for  one  part  of  the  time  range  than  another  as  deemed  advisable.  With  all  N 
particle  states  known  at  time  t^,  for  each  particle  in  turn  a  Monte  Carlo 
procedure  is  used  to  compute  the  further  history  up  to  time  t^p  +  ^ .  First, 
a  free-flight  time  t  is  chosen  for  the  particle  and  (t^p  +  t)  is  compared 
with  t^p  +  The  particle  trajectory  is  then  computed  up  to  the  lesser  of 
these  two  times.  If  a  scattering  occurs  within  the  time  interval  At,  the 
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result  of  the  col  lisle  and  the  new  flight  time  following  the  collision  are 
obtained  by  statistic;-  selection  procedures  in  the  program.  The  process  is 
repeated  until  t^p  +  '  is  reached.  It  is  unnecessary  to  continue  the  same 
free  flight  from  befor.  in  observation  time  to  after.  The  particle  must  have 
the  same  wave  vector  k  nd  positic  r  just  after  the  observation  time  t^  as 

just  before,  but  it  ct.:  oe  assigned  a  new  flight  time  commencing  at  t^. 

Thus,  the  remainder  of  le  flight  duration  at  the  observation  time  t^  (on 
the  flight  which  the  p.  tide  was  following  when  it  reached  t^)  need  not  be 
stored  for  each  partic'  .  The  values  of  the  particle  variables  at  time  t^ 
are  stored  for  all  N  -  ‘ides,  and  transport  quantities  or  interparticle 
effects  are  calculate-  ir  this  time.  Thus,  the  scatterings  and  flights  are 

computed  in  essential-,  .he  same  way  as  in  SPMC,  except  that  the  basic  segment 

of  a  particle's  histor;  is  the  interval  At  between  observation  times  rather 
than  the  (variable)  ir  -val  between  scatterings.  Each  segment  concludes  at 
an  observation  time,  ar  the  next  one  begins  at  this  time.  The  computer  time 
required  will  depend  m:  ily  on  the  number  of  scatterings  rather  than  the  num¬ 
ber  of  observation  tiir.-.  . 

9.5  ENERGY  BAND  STRUCT. TE 

For  a  cubic  type  ; _ni conductor,  the  momentum  vectors  at  the  energy  minima 
of  a  conduction  band  '  located  at  the  center  r  (p  *  0)  and  on  the  symmetry 
axes  (a  along  the  [1,  .  0]  directions  and  A  along  the  [I,  1,  1]  directions) 

of  the  first  Brillouir.  )ne.  Near  a  minimum  point  ])„,  the  E(p)  energy  function 
can  be  approximated  b}  >rms  to  second  order  of  a  Taylor  expansion  about  p0, 
i.e.,  by  a  quadratic  *  tion  of  ("p-"p0).  The  quadratic  form  is  called  the 
parabolic  approximate  o  Efp)  or  the  effective  mass  approximation.  The 

corresponding  surface:  constant  energy  centered  at  r  are  spheres  while 

those  centered  at  a  nr.  uim  on  a  a  or  A  axis  are  generally  ellipsoids  of 
revolution  having  a  cc  in  symmetry  axis  along  the  particular  a  or  A  axis. 

The  two  upper  valence  ds  have  degenerate  maxima  at  r.  The  equienergetic 
surfaces  about  r  are  ’  second  order  in  the  components  of  p,  more  complicated 
warped  surfaces. 

Further  up  the  e  ;y  valley  centered  at  the  minimum  point  p0  of  a 

conduction  band,  the  gy  function  deviates  from  a  simple  quadratic  expres¬ 

sion.  An  energy-momen  i  relation  accurate  to  larger  values  of  (p-p0),  into 
the  region  of  nonparat  city,  is  obtained  through  the  relationship 
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ECp)  =  E ( 1  +  ciE),  where  E(p)  on  the  left  denotes  the  quadratic  function  of 
momentum  characterizing  the  parabolic  approximation  and  a  is  the  nonparaboli- 
city  parameter.  This  device  accounts  partly  for  the  nonparabolicity  of  the 
band  structure  at  energies  further  up  the  valley. 

These  relatively  simple  approximations  of  band  structure  are  adequate 
for  electron  energies  up  to  about  0.5  eV  above  a  valley  minimum.  They  have 
been  widely  used  in  transport  calculations.  In  Monte  Carlo  simulation  they 
give  agreement  with  experiments  for  electric  fields  up  to  approximately 
10  kV/cm.  The  values  of  the  parameters  entering  into  the  expressions  (effec¬ 
tive  masses,  energy  values  at  extrema,  location  of  energy  minima  on  a  axes, 
nonparabolicity  parameters,  valence-band  warped  surface  parameters)  are  well 
established  for  all  the  important  cubic  semiconductors  of  Group  IV  elemental 
type  and  III-V  compound  type. 

The  present  form  of  the  ensemble  Monte  Carlo  program  allows  any  analytic 
approximation  for  the  energy  band  function  to  be  used.  This  has  been  accom¬ 
plished  with  relatively  simple  programming  compared  to  the  method  used  here¬ 
tofore  for  ellipsoidal  energy  surfaces  (Ref.  23).  The  algorithm  used  exploits 
the  high  symmetry  of  the  first  Brillouin  zone  in  cubic  crystals. 

The  set  of  momentum  vectors  specifying  the  Bloch  states  of  an  energy  band 
can  be  chosen  to  lie  in  the  volume  of  momentum  space  called  the  first 
Brillouin  zone.  The  zone  is  centered  at  the  origin  and  has  symmetry  proper¬ 
ties  reflecting  the  spatial  symmetry  of  the  crystal  lattice.  For  a  cubic 

crystal,  it  has  the  symmetry  of  a  cube.  To  any  momentum  vector  p  there 

correspond  47  other  momentum  vectors  in  the  Brillouin  zone  having  exactly 

the  same  energy  value  E(p).  The  48  vectors  are  related  by  the  48  symnetry 

operations  of  the  cube. 

This  synwetry  of  the  band  structure  can  greatly  facilitate  the  use  of 
analytic  expressions  corresponding  to  extremum  points  in  a  relatively  small 
part  of  the  Brillouin  zone,  namely,  in  an  irreducible  wedge  comprising  (1/48) 
of  the  total  volume.  To  determine  the  energy  value  for  an  arbitrary  momentum 
vector  of  the  Brillouin  zone,  it  is  transformed  by  the  appropriate  symmetry 
operation  into  the  corresponding  vector  of  the  irreducible  wedge.  The 
analytic  approximation  for  the  energy  about  the  nearest  extremum  point  in  the 
wedge  is  then  used  to  calculate  the  energy. 
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9.6  REALISTIC  BAND  STRUCTURE 

At  electron  energies  of  1  eV  or  more  above  the  minimum  point  of  an  energy 
valley  in  the  conduction  band,  the  approximate  analytic  expression  for  the 
energy  function  is  very  inaccurate.  Hence,  for  electric  fields  at  which 
average  electron  energies  exceed  0.5  eV,  more  accurate  band  structure  becomes 
necessary.  Furthermore,  at  even  lower  electron  energies,  interband  transitions 
to  the  second  conduction  band  begin  to  have  an  effect.  Then  the  Monte  Carlo 
simulation  must  include  interband  as  well  as  intraband  scattering  involving 
the  first  two  or  more  conduction  bands. 

For  electric  fields  above  100  kV/cm,  the  number  of  carriers  whose  energy 
exceeds  the  valence-conduction  bandgap  starts  to  be  significant.  Another 
carrier  collision  process,  impact  ionization,  then  becomes  important.  Impact 
ionization  involves  relatively  few  carriers  in  the  distribution,  those 
occupying  the  comparatively  small  high-energy  tail.  Calculation  of  ionization 
rates,  however,  requires  a  fairly  precise  determination  of  their  energy  dis¬ 
tribution.  Hence,  realistic  band  structure  is  a  necessity. 

Energy  bands  can  be  calculated  precisely  by  various  advanced  methods. 

Of  these,  the  empirical  pseudopotential  method  has  been  found  to  be  the  most 
useful  for  transport  problems.  Shichijo  and  Hess  (Ref.  23)  were  the  first  to 
employ  realistic  band  structure  in  Monte  Carlo  simulation.  Band  structure 
synmetry  also  facilitates  the  use  of  computed  energy  values,  since  the  energy 
and  momentum  gradient  of  energy  need  be  tabulated  only  at  a  mesh  of  points  in 
an  irreducible  wedge.  Again,  an  arbitrary  momentum  vector  in  the  Brillouin 
zone  is  transformed  into  the  corresponding  vector  in  the  wedge.  The  desired 
energy  value  is  obtained  from  a  weighted  average  of  the  energies  at  the  mesh 
points  surrounding  the  transformed  momentum  vector. 

9.7  SCATTERING  RATES 

In  addition  to  the  energy-band  structure,  Monte  Carlo  simulation  of 
charge  transport  in  a  semiconductor  crystal  requires  the  scattering  rates  of 
the  various  carrier  interaction  processes.  Both  properties  of  the  material 
must  be  known  fairly  precisely  at  all  carrier  energies  of  importance  in  the 
distribution  function  for  the  field  strength  considered.  In  general,  this 
becomes  more  difficult  as  field  strength  increases  and  produces  higher 
carrier  energies. 
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Many  high  field  effects  occur  at  electric  fields  F  <  50  kV/cm.  Within 
this  range  it  is  usually  adequate  to  treat  electron  transport  with  only  the 
first  conduction  band  and  approximate  analytic  energy  functions.  Scattering 
rates  are  calculated  using  ordinary  perturbation  theory  with  the  energy  density 
of  states  given  by  the  approximate  band  structure.  The  resulting  angular 
dependence  and  total  scattering  rates  of  the  various  scattering  mechanisms 
have  been  extensively  applied  in  Monte  Carlo  calculations. 

They  give  good  agreement  with  measured  high-field  transport  behavior  at 
electric  fields  below  50  kV/cm,  where  the  range  of  significant  electron 
energies  lies  below  1  eV.  Thus,  the  low  energy  scattering  rates,  valid  for 
electron  energies  below  1  eV,  are  now  well  established  for  the  important 
semiconductor  materials. 

Analytical  expressions  giving  the  scattering  rates  in  the  low  energy 
approximation  are  available  in  the  literature  for  all  the  important  scattering 
mechanisms.  These  include  intravalley  scattering  by  acoustic  phonons  and  by 
polar  and  nonpolar  optical  phonons.  In  intervalley  scattering,  the  electron 
is  scattered  from  one  valley  to  another  in  the  same  band,  either  an  equivalent 
valley  or  an  inequi valent  one  (the  latter  has  a  minimum  energy  different  from 
that  of  the  initial  valley).  Intervalley  phonon  scattering  is  induced  by 
acoustic  or  nonpolar-optical  phonons  of  rather  large  energy.  It  becomes 
important  at  electron  energies  high  enough  for  emission  of  these  phonons  to 
be  energetically  possible.  In  piezoelectric  crystals,  intravalley  scattering 
by  acoustic  piezoelectric  phonons  can  be  appreciable. 

Nonphonon  scattering  mechanisms  often  Important  are  ionized  impurity  and 
neutral  Impurity  scattering.  The  latter  can  also  cause  intervalley  scattering. 
In  ternary  alloys,^«Bi-«C  Intravalley  alloy  scattering  must  be  considered 
when  the  mole  fraction  x  of  the  binary  AC  is  not  close  to  the  value  zero  or 
one.  Alloy  scattering  can  be  viewed  as  arising  from  randomly  distributed 
potential  wells  against  a  uniform  periodic  background.  These  wells  occur  at 
the  lattice  sites  associated  with  one  of  the  allowing  constituents,  say  B, 
and  represent  the  difference  In  Ionic  potentials  at  sites  occupied  by  ion 
cores  of  A  and  B. 

In  Impact  Ionization  the  interaction  of  a  mobile  carrier  with  a  valence 
electron  causes  a  transition  of  the  latter  to  the  conduction  band,  so  the 
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final  state  has  an  electron-hole  pair  in  addition  to  the  primary  carrier.  The 
process  has  an  energy  threshold  that  is  not  less  than  the  bandgap  energy.  It 
is  important  for  fields  exceeding  100  kV/cm,  where  the  number  of  electrons 
with  energy  above  threshold  becomes  significant. 

A  scattering  function  calculated  by  perturbation  theory  is  proportional 
to  the  energy  density  of  final  states.  The  density  of  states  obtained  from 
the  effective  mass  approximation  is  customarily  inserted  into  the  calculation. 
However,  this  approximate  density  continually  increases  with  energy,  whereas 
the  true  density  of  states  peaks  at  some  energy  Ec  and  thereafter  falls  off. 

The  true  density  can  be  taken  as  that  given  by  the  realistic  band  structure 
calculation.  For  electron  energies  below  E  ,  the  true  and  approximate  densi- 

v 

ties  increase  with  energy  and  agree  over  much  of  the  interval.  The  perturbation- 
theoretic  scattering  rates  of  nonpolar  phonon  processes  as  well  as  the  total 
scattering  rate  increase  with  energy  below  E  in  agreement  with  established 
values  for  the  rates.  However,  these  rates  continue  to  increase  indefinitely 
with  energy,  contrary  to  more  accurate  calculations  which  show  a  fall  off  in 
the  vicinity  of  Ec.  In  practice  the  scattering  rates  from  perturbation  theory, 
with  coupling  parameters  adjusted  to  give  agreement  with  observed  transport 
behavior  up  to  the  lower  range  of  high  fields,  can  be  used  for  electron  ener¬ 
gies  below  E_.  At  higher  energies  the  scattering  rates  are  modified  in 

w 

accordance  with  calculated  rates  that  are  based  on  non-perturbatlve  field- 
theoretic  methods  and  use  the  true  density  of  states  (Ref.  23  and  24). 

9.8  QUANTUM  EFFECTS 

At  electron  energies  about  1.5  eV  above  the  conduction  band  edge,  quantum- 
mechanical  effects  become  Important  in  semiconductor  transport.  At  these 
energies,  the  number  of  electron  states  per  unit  energy  range  is  high  and  the 
scattering  rates  of  many  phonon  collision  processes  are  large.  As  dictated 
by  the  energy-time  uncertainty  effect,  scattering  introduces  uncertainty  into 
the  electron  energy.  At  the  highest  scattering  rates  the  energy  uncertainty 
is  not  necessarily  small  compared  to  the  energy  itself. 

In  the  quasi-classical  Boltzmann  equation  description,  the  energy  eigen¬ 
states  of  the  unperturbed  electron  are  used  in  defining  the  distribution 
function.  But  when  the  electron  interacts  strongly  with  lattice  phonons,  it 
is  more  accurately  described  as  propagating  in  a  perturbed  state.  In  a 
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perturbed  state  the  electron  continually  emits  phonons  and  reabsorbs  them  by 
virtual  collision  processes.  As  a  result  of  the  interactions,  the  state  has 
a  range,  aE,  of  energy  values  and  an  energy  center  displaced  from  the  sharp 
energy  value  of  the  corresponding  Bloch  state. 

Much  of  the  incorrect  nature  of  the  quasi -classical  Boltzmann  equation 
can  be  removed  by  introducing  the  concepts  of  a  quasi-particle  and  its  self¬ 
energy.  The  quasi-particle  is  the  Bloch  electron  with  its  accompanying  cloud 
of  virtual  phonons.  The  self-energy  corresponds  to  a  mass  renormalization  and 
can  be  calculated  by  field-theoretic  operator  techniques.  The  real  part  of 
the  self-energy  is  equal  to  the  level  shift  and  the  imaginary  part  is  equal  to 
the  collision-broadened  width,  aE,  of  the  level.  Also,  the  total  scattering 
rate  (in  real  collisions)  of  the  quasi-particle  is  directly  proportional  to  its 
level  width,  aE,  from  collision  broadening,  in  accordance  with  the  uncertainty 
effect. 

Recently  it  has  been  shown  by  H.  Shichijo  and  K.  Hess  (Ref.  23)  and  by 
J.Y.  Tang  and  K.  Hess  (Ref.  24)  for  gallium  arsenide  and  silicon,  respectively, 
that  the  accuracy  of  impact  ionization  rates  calculated  by  the  Monte  Carlo 
method  Is  significantly  improved  by  introducing  the  quantum  effects  described 
above.  In  the  simulation,  the  particle  is  scattered  into  a  range  of  final 
energies  corresponding  to  the  widths  of  the  initial  and  final  quasiparticle 
states.  Moreover,  the  total  scattering  rate  at  energies  above  Ec  is  chosen 
in  accordance  with  the  energy  dependence  of  the  rate  given  by  the  self-energy 
calculation.  For  carrier  energies  below  Ec  the  perturbation-theoretic 
scattering  rates  are  used,  and  the  high  energy  rates  are  joined  continuously 
onto  these. 

9.9  TIME-OF-FLIGHT  ALGORITHMS 

A  free-flight  time  is  the  interval  between  successive  collisions,  during 
which  the  electron  drifts  under  the  influence  of  the  field.  The  algorithm 
that  selects  the  flight  times  is  one  of  the  most  important  features  of  a  Monte 
Carlo  computer  program.  A  better  method  improves  the  speed  of  the  program  and 
the  accuracy  of  the  simulation.  For  our  purposes,  the  procedure  must  be 
relatively  easy  to  implement  even  with  a  complicated  energy-band  stucture, 
such  as  that  seen  by  the  higher  energy  electrons  of  the  distribution. 
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Another  consideration  is  compatibility  with  spatial  dependence.  In  the 
operation  of  many  devices,  carrier  transport  varies  with  position  in  space. 
This  arises  from  material  nonuniformity  and  interfaces  or  from  prescribed 
boundary  conditions,  such  as  carrier  injection  at  a  junction  or  a  hetero¬ 
junction  barrier.  Often  a  nonuniform  field  occurs.  Hence,  it  is  desirable 
that  the  algorithm  be  adaptable  to  the  case  of  a  rapidly  varying  field. 

Let  v{t)  be  the  total  collision  rate  of  an  electron  at  time  t  along  the 
path  of  a  free  flight.  The  probability  of  the  electron  being  scattered  in  the 
time  interval  (M  +  A«)  can  be  taken  as 

P»c (t,  St)  =  i/(t)  At  (46) 

provided  At  is  short  enough  for  W<)  to  be  essentially  constant  over  the 
interval.  It  follows  directly  that  the  probability  of  the  electron  drifting 
without  collision  from  t=0tot=T  is 

P{T)  =  exp  -  f  (i/p{t))dt  ( 47 ) 

The  cumulative  probability  that  the  free-flight  time  lies  somewhere  in  the 
time  interval  (0,  T)  is  1  -  P(T),  which  has  the  value  0  at  T  =  0  and 
approaches  the  value  1  for  large  T.  Free-flight  times  having  the  required 
probability  distribution  can  be  selected  as  follows:  Generate  a  random  number 
r,  uniformly  distributed  between  0  and  1  and  solve 

r  *  1  -  P(T)  (48) 

for  t,  which  is  taken  to  be  the  free-flight  time.  However,  the  collision  rate 
is  a  complicated  function  of  the  time  along  the  electron's  free  trajectory, 
and  there  is  no  direct  way  of  inverting  Equation  48  to  obtain  T.  Each  selec¬ 
tion  of  a  flight  time  involves  a  numerical  integration,  as  shown  in  Equation 
47. 

Two  alternative  methods  are  available.  The  most  commonly  used  device  is 
to  introduce  self-scattering,  a  fictitious  process  which  leaves  the  electron 
state  unchanged  but  serves  to  make  the  total  scattering  rate  (including  self¬ 


scattering)  a  constant  for  all  electron  states,  p.  Then,  the  right-hand  side 
of  Equation  47  becomes  a  simple  exponential  function  of  T  and  Equation  48, 
and  is  easily  inverted  to  express  T  as  a  logarithmic  function  or  r.  The 
drawback  to  the  self-scattering  scheme  is  evident  at  energies  where  the  total 
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of  the  real  processes  is  relatively  small  and  the  self-scattering  rate  is  then 
necessarily  large.  As  a  result,  an  inordinate  amount  of  computer  time  is  spent 
on  self-scattering  events  which  have  no  physical  significance.  To  improve  the 
efficiency  of  the  simulation,  one  can  divide  momentum  space  into  energy  regions 
and  assign  optimum  values  to  the  total  scattering  rate  in  each  region.  The 
method  works  well  if  simple  analytical  expressions  for  the  energy  surfaces 
can  be  used,  as  is  the  case  at  field  strengths  where  the  electrons  are  concen¬ 
trated  mainly  in  the  lower  regions  of  the  energy  valleys.  But  when  realistic 
band  structure  is  required,  as  at  high  electron  energies,  the  subdivision  into 
energy  regions  becomes  intractable. 

The  most  utilitarian  device  dispenses  with  the  integral  form  of  Equation 
48  and  utilizes  the  differential  form  of  Equation  46.  The  procedure  is  the 
following.  Consider  an  electron  in  free  flight  at  time  t.  First  calculate  the 
total  scattering  rate  t/( t )  for  the  state  of  the  electron  at  this  time.  The 
random  number  4,  uniformly  distributed  between  0  and  1,  is  generated  to  decide 
whether  the  electron  is  scattered  in  the  next  time  increment  At.  If  r  satis¬ 
fies  r  <  i/(t)At  proceed  to  scatter  the  electron  at  t  +  [r/i/(t)].  Otherwise, 
continue  its  free  flight  to  time  Wt,  at  which  time  the  entire  process  is 
repeated.  Experience  reported  with  this  method  (Ref.  24)  indicates  that  as 
long  as  At  is  chosen  to  be  about  one  tenth  the  mean  scattering  time  in  the 
current  region  of  momentum  space,  both  the  speed  and  the  accuracy  of  the 
computation  are  satisfactory.  The  differential  scheme  has  the  following 
advantages:  (I)  no  complications  arise  from  the  complexity  of  the  energy 
surfaces,  and  (2)  a  simple  generalization  makes  it  applicable  to  spatially 
varying  fields  (where  the  integral  method  is  unsuitable). 

9.10  COLLECTIVE  EFFECTS 

When  carrier  concentrations  are  high,  charge  transport  can  be  strongly 
influenced  by  collective  effects.  These  arise  from  interaction  among  the 
particles  and  the  resulting  many-particle  correlations.  It  is  convenient  to 
resolve  the  interaction  between  carriers  into  two  parts.  One  part,  described 
as  carrier-carrier  scattering,  accounts  for  the  strong  short-range  coulomb 
interaction  between  an  individual  pair  of  particles  in  the  immediate  vicinity 
of  one  another.  Each  carrier-carrier  collision  conserves  the  total  energy  of 
the  pair  but  exchanges  energy  between  them.  The  collisions  are  effective  in 
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equilibrating  the  carrier  energy  distribution.  This  type  of  scattering  must 
be  included  in  the  Monte  Carlo  simulation  at  large  carrier  densities  to  des¬ 
cribe  the  strong  tendency  toward  thermal izati on  of  the  carrier  distribution 
under  these  conditions.  Clearly  a  many-particle  ensemble  Monte  Carlo  scheme 
is  required. 

The  other  part  of  the  carrier-carrier  interaction  is  included  in  the 
space-charge  effect  and  describes  the  average  continuous  influence  on  any 
single  carrier  of  all  the  others.  This  accounts  for  the  long-range  part  of 
the  coulomb  interaction  and  contributes  to  the  total  field  acting  on  a  carrier. 
The  field  in  the  material  is  the  result  of  both  the  external  applied  potentials 
and  the  spatial  distribution  of  charge  from  both  the  carriers  and  the  back¬ 
ground  of  ionized  impurities.  The  resulting  self-consistent  electric  field 
can  be  calculated  by  Poisson's  equation  when  the  charge  distribution  is  known. 

The  MPEMC  scheme  is  well  suited  to  the  analysis  of  collective  effects. 

As  an  example,  consider  carrier-carrier  collision  (Ref.  22).  Two  carriers 
enter  a  mutual  scattering  in  the  initial  states  (p^,  "pO  and  emerge  in  the 
final  states  (Pj',  p2')  with  probability  equal  to  S2Cp^,  P2;  P^»  P2').  where 
S2  is  the  pair  scattering  function.  The  final  states  ("Pj%  ”p2')  are  those 
allowed  by  energy  and  crystal -momentum  conservation.  Integrating  S2  over  the 
entire  fHspace  of  the  variables  Pj',  p^  gives  t2  (j^,  p2),  the  total  mutual 
scattering  rate  for  the  initial  pair  ("p ^ ,  ~p2 ) .  For  carrier  1,  the  scattering 
probability  per  unit  time  for  this  process  at  time  t  is  proportional  to 
T2  (plt  ”p2)  f  (p2,  r,  t).  Each  carrier  (and  each  impurity  ion)  has  its 
coulomb  field  screened  by  local  adjustment  of  the  mobile  charges  in  its 
vicinity.  Hence,  S2  describes  scattering  by  a  screened  coulomb  potential 
(i.e.,  similar  to  charged  impurity  scattering).  The  microscopic  details  of 
the  screening  depend  on  the  carrier  distribution  function  in  the  neighborhood 
of  the  charge  and  an  exact  treatment  is  complicated  in  nonequilibrium 
situations. 

The  carrier-carrier  scatterings  are  calculated  using  the  distribution 
function  obtained  at  the  preceding  observation  time.  In  a  situation  of  spatial 
uniformity,  the  momentum  distribution  of  neighboring  carriers  can  be  estimated 
by  sampling  from  the  states  of  all  the  other  carriers  in  the  ensemble  at  the 
observation  time.  The  observation  times  may  be  made  as  frequent  as  necessary 
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for  this  purpose  with  an  addition  to  computer  time  that  does  not  include  addi¬ 
tional  single-particle  scattering  computations.  The  screening  is  treated 
approximately  by  Debye  theory,  which  is  valid  in  equilibrium.  Simulation  of 
carrier-carrier  scatterings  over  the  segment  of  each  particle  history  that 
occurs  in  the  time  interval  between  t^  and  t^m  +  ^  can  be  done  with  the 
following  successive  steps: 

1.  Consider  particle  1  in  free  flight  in  this  interval.  Choose  a 

companion  state  "p.  randomly  from  the  other  N-l  particle  states  in  memory 
J  /  —  \ 

storage  at  the  observation  time  V  . 

2.  Consider  any  statistical  selection  involved  in  choosing  the  dura¬ 
tion  or  the  type  of  terminating  collision  for  this  flight.  Use  the  pair 
(?1»  Pj )  as  the  initial  states  and  r2  (p^  p^)  as  the  scattering  probability 
of  a  mutual  scattering  by  particle  1,  where  p^  is  the  momentum  of  particle  1 
at  the  time  in  the  free  flight  at  which  the  statistical  selection  is  being 
made.  (See  Section  9.9.)  Of  course,  the  competing  single-particle  scatter¬ 
ing  probabilities  used  at  this  time  are  those  evaluated  at  "p, . 

3.  If  Step  (2)  results  in  a  mutual  scattering,  then  select  values  of 
(p*: ",  "pj")  in  accordance  with  the  distribution  given  by  S2  (Pj,  Pj ;  p^',  Pj')- 

4.  Repeat  Steps  (1)  to  (3)  for  all  the  flight/scattering  cycles 
required  for  particle  1  to  reach  the  observation  time  t^m  +  If  the 
probability  of  a  mutual  scattering  is  relatively  small  compared  to  the  other 
processes.  Steps  (1)  to  (3)  may  be  executed  in  only  a  fraction  of  the  flight/ 
scattering  cycles  while  replacing  T2  in  Step  (2)  by  T2/a. 

5.  Repeat  Steps  (1)  to  (4)  for  particles  2,  3,  ....  N. 

Due  to  the  long  range  of  the  coulomb  force,  each  charge  carrier  feels, 
to  some  extent,  the  presence  of  all  other  carriers  and  all  impurity  ions  in 
the  device.  These  other  charges  contribute  to  the  total  field  (beyond  that 
of  the  perfect  ionic  lattice)  acting  on  the  carrier  in  free  flight.  Because 
of  the  large  number  involved,  their  net  contribution  to  the  field  acts  in  a 
continuous  manner  and  is  determined  by  the  macroscopic  spatial  distribution  of 
charge,  the  sum  of  the  charge  distributions  given  by  the  carrier  density 
function  and  the  density  distribution  of  impurity  ions. 


The  space  charge  effect  can  be  included  in  the  computation  when  particle 
positions  are  recorded.  The  carrier  density  function  at  a  particular  time  is 
estimated  from  the  location  of  particles  in  the  carrier  ensemble  at  that  time. 
The  total  self-consistent  potential  at  t^  is  calculated  by  solving  Poisson's 
equation  with  the  space  charge  distribution  at  t^  and  the  prescribed 
potentials  applied  at  the  spatial  boundaries.  In  the  time  step  to  the  next 
observation  time  t^m  +  all  free  flights  are  presumed  to  be  governed  by 
this  potential.  At  the  next  observation  time,  the  self-consistent  field  is 
recalculated. 

9.11  NEGATIVE  DIFFERENTIAL  CONDUCTIVITY 

Differential  conductivity  is  the  rate  of  change  of  current  density,  J, 
with  electric  field,  F.  Negative  differential  conductivity  (NDC)  can  arise 
from  various  mechanisms  in  the  bulk  or  in  connection  with  junction  phenomena. 
It  shows  two  distinct  forms  of  behavior. 

In  voltage-controlled  manifestations  of  NDC,  the  curve  of  J  versus  F  has 
a  shape  resembling  an  N.  Beyond  the  threshold  value  of  F,  the  curve  has  a 
negative  slope  and  NDC  continues  to  some  larger  value  of  F  beyond  which  the 
slope  turns  positive  again.  Bulk  NDC  of  this  type  occurs  at  high  fields  in 
GaAs  and  other  polar  semiconductors  as  a  transferred-electron  effect.  With 
a  short  length  of  material  in  the  appropriate  external  circuit,  a  transient 
instability  occurs  in  which  a  region  of  high  field  and  charge  accumulation 
forms.  The  region  has  the  structure  of  a  layer  transverse  to  the  electric 
field  and  propagates  at  the  speed  of  the  applicable  drift  velocity.  As  soon 
as  the  layer  disappears  at  the  anode  side,  a  new  layer  forms  at  the  cathode 
side.  The  domain  propagation  is  utilized  in  the  Gunn  oscillator. 

In  current-controlled  NDC,  the  curve  of  J  versus  F  has  the  approximate 
shape  of  an  S  and  J  is  a  multivalued  function  of  F.  The  differential  con¬ 
ductivity  at  some  point  starts  to  rise  with  increasing  field.  At  some  break¬ 
over  value  of  F,  an  instability  in  J  develops  and  a  small  further  increase 
in  J  is  accompanied  by  a  switch  from  a  high-  to  a  low-field  region.  This  form 
of  NDC  leads  to  the  formation  of  a  high-current  filament  running  along  the 
field  direction.  The  constriction  of  the  current  results  in  local  heating  of 
the  material  and  further  increase  in  current.  If  the  current  is  allowed  to 
increase  in  this  breakdown  regime,  melting  in  the  filament  region  will  damage 
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the  device.  Analysis  has  begun  using  Monte  Carlo  simulation  to  determine  if 
current-controlled  NDC  can  occur  in  the  bulk  from  a  transferred  electron 
effect.  This  could  provide  basic  understanding  of  second  breakdown  (Refs.  1, 
25  and  26),  the  damage  mechanism  commonly  observed  at  very  high  fields. 


10.0  MONTE  CARLO  ANALYSES  AT  VERY  HIGH  FIELDS 


Monte  Carlo  simulation  of  nonlinear  transport  at  fields  in  the  range  from 
100  to  1000  kV/cm  have  been  performed  in  collaboration  with  Professor  K.  Hess 
and  T.  Wang  at  the  University  of  Ill inois-Urbana,  Illinois,  under  subcontract 
to  Rocketdyne*.  The  analysis  addresses  some  basic  questions  concerning  the 
physics  of  second  breakdown  in  bulk  semiconductors  and  is  the  first  accurate 
modeling  of  its  transport  behavior.  Professor  Hess  has  pioneered  in  the 
extention  of  Monte  Carlo  calculations  to  high  fields.  His  group  has  developed 
and  applied  computer  programs  containing  many  of  the  refinements  required, 
particularly  realistic  band  structure  and  accurate  scattering  rates  at  high 
energies.  Hence,  his  interest  and  involvement  in  the  study  has  enabled  us 
to  carry  out  reliable  calculations  within  the  time  limits  set. 

The  numerical  simulations  were  performed  for  electron  transport  in  GaAs, 
Si  and  InP  and  used  the  most  complete  Monte  Carlo  program  presently  in 
existence.  The  calculations  simulate  electron  motion  and  multiplication 
(impact  ionization),  the  empirical  pseudopotential  method.  All  scattering 
mechanisms  relevant  at  high  energies  are  included  in  the  simulation.  Band- 
structure  changes  with  temperature  are  also  taken  into  account.  The  basic 
questions,  raised  by  M.E.  Snyder,  are  the  following: 

1.  Can  electrons  escape  to  higher  conduction  bands  exhibiting  higher 
mobility?  In  other  words,  is  it  possible  that  at  very  high  electric  fields 
and/or  during  short  items  the  current  does  not  saturate  or  decrease  but 
increases  with  electric  field  even  if  the  carrier  concentration  were  to 
remain  constant? 

2.  Is  it  possible  for  high  level  current  injection  to  be  augmented  by 
ballistic  transport  during  short  times? 

3.  Does  avalanche  breakdown,  as  often  conjectured  previously,  generally 
die  out  as  the  temperature  increases? 

4.  Is  the  temperature  dependence  of  negative  differential  conductivity 
(NDC)  adequate  enough  to  support  microplasma  formation  and  growth  alone? 


*  Rockwell  International  Corporation 
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10.1  RESULTS  FOR  ITEM  1 

The  answer  is  yes.  A  clear  increase  in  drift  velocity  can  be  seen  for 
GaAs  (Figs.  A4,  A7)  and  for  InP  (Figs.  A13,  A16)  at  higher  temperatures.  But 
no  increase  is  computed  for  Silicon  (Fig.  Al)  up  to  10,000  kV/cm.  Thus,  the 
process  appears  to  be  very  complicated  and  highly  dependent  on  the  material, 
i.e.,  on  specific  band  structure  and  scattering-rate  details. 

10.2  RESULTS  FOR  ITEM  2 


Typical  scattering  rates  at  high  energies  as  given  by  the  numerical 
simulations  were  investigated.  The  results  show  that  a  possibility  exists 
for  ballistic  motion  after  the  electrons  have  attained  sufficient  energy  to 
overcome  the  maximum  in  the  density  of  states  (DOS),  (Fig.  A2).  At  these 
high  energies  ballistic  motion  may  persist  for  picoseconds,  while  at  the 
energy  of  highest  density  of  states  the  ballistic  transport  is  limited  to 
less  than  0.01  ps. 

10.3  RESULTS  FOR  ITEM  3 

Surprisingly,  the  result  is  no.  Avalanche  breakdown  does  not  die  out 
as  the  temperature  increases.  It  only  becomes  weaker  in  a  small  temperature 
range  up  to  about  400  K  or  500  K  in  the  case  of  GaAs  (Fig.  A3).  At  higher 
temperatures  the  rate  becomes  stronger  (Figs.  A5,  A12,  715).  The  reason  seems 
to  be  a  quantum  mechanical  renormalization  which  does  not  let  the  scattering 
rates  grow  with  temperature  much  above  their  peak  value  at  the  point  of  the 
highest  DOS.  It  is  the  almost  constant  peak  value  which  is  responsible  for 
the  ionization  rate.  This  result,  however,  is  preliminary  and  needs  to  be 
carefully  checked. 

10.4  RESULTS  FOR  ITEM  4 

The  answer  is  no,  as  can  be  seen  from  the  magnitude  of  the  temperature 
variations  in  Sections  10.1  and  10.3. 


•4 


66 


AFWL-TR-85-116 


REFERENCES 


1.  Snyder,  M.E.,  The  Physics  of  Second  Breakdown,  AFWL-TR-83-22 ,  Air  Force 
Weapons  Laboratory,  Kirtland  Air  Force  Base,  New  Mexico,  pp.  7-8,  Nov  1983. 

2.  Bower,  R.,  and  Steele,  M.C.,  "Plasma  Effects  in  Solids,"  Proc.  of  the 
IEEE,  Vol .  52,  p.  1105,  Oct  1964. 

3.  Platzman,  P.M.,  and  Wolff,  P.A.,  "Waves  and  Interactions  in  Solid  State 
Plasmas,"  Solid  State  Physics,  Supp.  13,  Academic  Press,  New  York,  New 
York,  1973. 

4.  Bloemberger,  N.,  et  al ,  "Laser- Induced  Electric  Breakdown  in  Solids," 

IEEE  Joul.  of  Quant.  Elec.,  Vol  QE-10,  No.  3,  p.  375,  Mar  1974. 

5.  Chen,  F.F.,  Introduction  to  Plasma  Physics,  Plenum  Press,  New  York, 

New  York,  pp.  1-4,  1974. 

6.  Bube,  R.H.,  Electronic  Properties  of  Crystalline  Solids,  Academic  Press, 
New  York,  New  York,  pp.  121-160,  1974. 

7.  Kittel,  C.,  Introduction  to  Solid  State  Physics,  Fifth  Edition,  John 
Wiley  and  Sons,  New  York,  New  York,  Table  2,  p.  223,  1976. 

8.  Bube,  R.H.,  IBID,  pp.  468-469. 

9.  Bennett,  W.H.,  "Magnetically  Self-Focusing  Streams,"  Phys.  Rev., 

Vol.  45,  p.  890,  Jun  1934. 

10.  Willardson,  R.K.,  et  al ,  "Plasmas  in  Semiconductors  and  Semimetals," 
Semiconductors  and  Semimetals,  Vol.  1,  Academic  Press,  New  York,  New 
York,  pp.  444-447,  1966. 

11.  Seshadri ,  S.R.,  Fundamentals  of  Plasma  Physics,  American  Elsevier 
Publishing  Company,  pp.  12$-i2/r,  1973. 

12.  Boyd,  G.D.,  et  al ,  "Excitation  of  Plasma  Oscillations  and  Growing  Plasma 
Waves,"  Phys.  Rev.,  Vol.  109,  p.  1398,  1958. 

13.  Pines,  D.,  and  Schrieffer,  J.R.,  "Collective  Behavior  in  Solid  State 
Plasmas,"  Phys.  Rev.,  Vol.  124,  No.  1,  p.  1387,  Dec  1961. 

14.  Willardson,  R.K. ,  et  al ,  IBID,  p.  478. 

15.  Ancker-Johnson,  B.,  and  Drunmond,  J.E.,  "Thermal  Pinching  in  Electron-Hole 
Plasma,"  Phys.  Rev.,  Vol.  131,  No.  5,  p.  1961,  Sep  1963. 

16.  Chen,  W.S.,  and  Ancker-Johnson,  B.,  "Pinch  Oscillations  In  Electron-Hole 
Plasmas,  I.  Theory  and  II.  Experiment,"  Phys.  Rev.  B.,  Vol.  2,  No.  11, 
pp.  4468  and  4477,  Dec  1970. 


67 


AFWL-TR-85-116 


REFERENCES  (Concluded) 


17.  A  general  survey  of  recent  developments  is  given  in  "The  Physics  of 
Submicron  Structures,"  edited  by  Grubin,  H.L.,  Hess,  K.,  Iafrate,  G.J. 
and  Ferry,  D.K.,  Plenum  Press,  New  York,  New  York,  1984. 

18.  Hess,  K.,  Advances  in  Electronics  and  Electron  Physics,  Vol .  59, 

pp.  239-291,  1982.  “ 

19.  Sze,  S.E.,  "Physics  of  Semi-conductor  Devices,"  Second  Edition,  Wiley- 
Interscience,  New  York,  New  York,  p.  298,  1981. 

20.  Grubin,  H.L.,  Ferry,  D.K.,  Iafrate,  G.J.  and  Barker,  J.R.,  VLSI 
Electronics:  Microstructure  Science,  Vol.  3,  Academic  Press,  New  York, 

New  York,  1982. 

21.  Jacoboni,  C.  and  Reggiani.  L..  Rev.  Mod.  Phys,  Vol.  55,  pp.  643-705, 

Jul  1983. 

22.  Price,  P.J.,  Semiconductors  and  Semimetals,  Vol.  14,  Academic  Press, 

New  York,  pp.  249-308,  1979. 

23.  Shichijo,  H.  and  Hess,  K. ,  Phys.  Rev.  B,  Vol.  23,  pp.  4197-4207,  Apr  1981. 

24.  Tang,  J.Y.  and  Hess,  K.J.,  Appl .  Phys.,  Vol.  54,  pp.  5139-5144,  Sep  1983. 

25.  Portnoy,  W.M.,  "Reverse  Bias  Second  Breakdown  in  Power  Switching 
Transistors,"  AFWL-TR-82-139,  Air  Force  Weapons  Laboratory,  Kirtland 
Air  Force  Base,  New  Mexico,  May  1983. 

26.  McMullen,  S.A.  and  Portnoy,  W.M.,  "Energy  Considerations  In  Second 
Breakdown,"  AFWL-TR-84-111,  Air  Force  Weapons  Laboratory,  Kirtland  Air 
Force  Base,  New  Mexico,  Feb  1985. 


68 


AFWL-TR-85- 1 16 


'1 


1 1* 
*»• 


APPENOIX 

GRAPHS  AT  UNIVERSITY  OF  IUINOIS-URBANA  MONTE  CARLO  RUNS 
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Figure  A2.  Density  of  states  versus  energy  in  silicons  conduction  band. 
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Figure  A5.  Average  electron  energy  versus  applied  field  in  FaAs  at  T  -  500  K. 
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Figure  A10.  Average  electron  velocity  versus  applied  field  for  InP  at  T 


Figure  All.  Average  electron  energy  versus  applied  field  in  InP  at  T 
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ng  constant  is  unchanged. 
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age  electron  velocity  versus  applied  field  in  InP  at  T  -  500  K. 
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electron  energy  versus  applied  field  in  InP  at  T  -  500  K. 
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gure  A16.  Average  electron  velocity  versus  the  applied  field  in  InP  at  T  =  800  K. 


AVERAGE  ELECTRON  ENERGY 


Average  electron  energy  versus  applied  field  in  InP  at  T  =  800  K. 


